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ABSTRACT 


The phenomenon known as aeroelastic divergence is the focus of this work. The 
analyses and experiment presented here show that divergence can occur without a 
structural dynamic mode losing its oscillatory nature. Aeroelastic divergence occurs 
when the structural restorative capability or stiffness of a structure is overwhelmed by the 
static aerodynamic moment. This static aeroelastic coupling does not require the 
structural dynamic system behavior to cease, however. Aeroelastic changes in the 
dynamic mode behavior are governed not only by the stiffness, but by damping and 
inertial properties. The work presented here supports these fundamental assertions by 
examining a simple system: a typical section airfoil with only a rotational structural 
degree of freedom. 

Aeroelastic stability analysis is performed in the discrete time domain. The aerodynamic, 
structural dynamic, and downwash relationships are cast as time-marching equations and 
combined to form aeroelastic state space equations. The discrete time eigenvalues and 
eigenvectors of the coupled system are computed. This method is advantageous because 
the exact roots and the degree of stability of the system are determined, within the 
framework of the aerodynamic and structural dynamic representations. The discrete-time 
eigenvalues are transformed into the continuous time domain to facilitate their 
interpretation. 

Results from the analysis have identified configurations of a simple model that exhibit 
different types of dynamic mode behavior as the system encounters divergence. For the 
simple configuration examined, these results indicate that low inertial properties and 
elastic axis location near the center of pressure promote divergence while the dynamic 
mode persists. Large inertias and large separation between elastic axis and center of 
pressure promote divergence where the dynamic mode becomes a static mode. 

A wind tunnel model was designed and tested to examine divergence experimentally. 

The experimental results validate the analytical calculations and explicitly examine the 
divergence phenomenon where the dynamic mode persists. Three configurations of the 
wind tunnel model were tested. The experimental results agree very well with the 
analytical predictions of subcritical characteristics, divergence velocity and behavior of 
the noncritical dynamic mode at divergence. 
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NOMENCLATURE 


Symbols: 

A Aerodynamic matrix 

a Amplitude of cycles, used in logarithmic decrement method 

B Aerodynamic matrix 

b Semi-chord 

C Aeroelastic force matrix 

Cux Lift curve slope 

D Structural dynamic matrix 

E Downwash matrix 

e Elastic axis position, measured positive aft from the center of 

pressure 

f Aerodynamic load vector 

f (x Frequency of torsional mode (Hz) 

I a Torsional inertia 

j Imaginary number, j = V~d 

k Reduced frequency, (k=a>b/U) 

Lwake Length of the wake 

M Number of aerodynamic elements on the wing 
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Mass 


N Total number of aerodynamic elements 

n Time step number 

Ncycles Number of cycles 

Nwake Number of aerodynamic elements in the wake 
q Generalized structural coordinates 

q Dynamic pressure 

r w Radius of gyration 

U Velocity 

V Reduced velocity, (V=U/(O a b) 
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x Vector of locations of vortices in the aerodynamic model; chord- 

wise location 

z Discrete time eigenvalue 

At Time step size, temporal discretization 

Ax Aerodynamic element size, spatial discretization 

T Vorticity vector 
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K (I Torsional stiffness 

M Moment 

a Angle of attack 
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5 Logarithmic decrement 

X Continuous time eigenvalue 

A Intermediate Southwell parameter, slope of moment vs. rigid angle 

of attack 
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p Density of air 
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n time step number 
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a pertaining to the torsional degree of freedom 

mode pertaining to a specified mode 
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jth column 
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wake quantity in the wake 

TE trailing edge 

A pertaining to the aerodynamics 

S pertaining to the structure 

0 vector or quantity at time 0 

1 matrix multiplies vector at time step n 

2 matrix multiplies vector at time step n+1 

o the rigid value of a quantity 
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INTRODUCTION 


Aeroelasticity is concerned with systems in which there is substantial interaction among 
the aerodynamic, inertial and structural forces of an object. The phenomenon known as 
aeroelastic divergence occurs when the structural restorative capability or stiffness of a 
structure is overwhelmed by the static aerodynamic moment. The static aeroelastic 
coupling that produces divergence does not require the dynamic system behavior to 
cease, however. Aeroelastic changes in the dynamic mode behavior are governed not 
only by the stiffness, but by damping and inertial properties. 

The work presented supports these assertions by examining a simple system: a typical 
section airfoil with only a rotational structural degree of freedom. The analyses and 
experiment to be presented show that divergence can occur without a structural dynamic 
mode losing its oscillatory nature and becoming static. 

The aeroelastic analysis method utilized in this study allows calculation of the 
eigenvalues or modal characteristics of a system for subcritical, critical and supercritical 
systems. The primary analysis is performed in the discrete time domain. The 
aerodynamic, structural dynamic, and downwash relationships are cast as time-marching 
equations and combined to form aeroelastic state space equations. The discrete time 
eigenvalues and eigenvectors of the coupled system are computed. The discrete-time 
eigenvalues are transformed into the continuous time domain to ease interpretation. This 
method is advantageous, as the exact roots and the degree of stability of the system are 
determined, to the extent of the accuracy of the aerodynamic and structural dynamic 
representations. The method differs from traditional aeroelastic analyses. Background 
information is provided on these traditional methods, which reveals differences between 
the current method and each of them. Most of the traditional analyses produce results 
that are only valid for neutrally stable behavior. This limitation is often not in the 
stability method itself, but rather in approximations in the aerodynamic behavior. The 
current analysis method resembles the p-method to be discussed, but the fundamental 
quantity, p, has been replaced by the discrete time unit delay operator, z. 

The discrete-time aeroelastic eigenanalysis method was used to examine the aerodynamic 
and structural parametric space of a typical section airfoil that had a single structural 
degree of freedom. These results distinguished configurations where different types of 
dynamic mode behavior were observed as the system encountered divergence. This 
facilitated the design of an experiment which encountered divergence while the structural 
dynamic mode persists. 

A wind tunnel model was designed and tested to examine divergence experimentally and 
validate the analytical calculations. All freedom of motion was denied to the airfoil. 



except for rotation about the elastic axis. Allowing only the single structural degree of 
freedom eliminated the complications of interpreting modal interaction effects or 
participation of multiple modes in the divergence mechanism. This simplicity allowed 
the focus to be precisely on the coupling of the aerodynamics with the structural pitching 
motion. Three configurations of the wind tunnel model were tested to examine the 
effects of a limited range of torsional stiffness and inertia. All three configurations 
exhibited divergence of a static mode existing simultaneously with a dynamic mode. The 
constraints imposed by the fundamental design of the model limit the potential source of 
both the statically unstable mode and the measured dynamic mode. The mode that 
originates as the structural dynamic pitch or torsional mode was tracked at subcritical 
airspeeds. As the dynamic pressure was increased, the aeroelastic coupling changes the 
damping and frequency of this tracked mode. At divergence, this mode appears as a 
damped oscillatory mode. The frequency and damping of the dynamic mode are 
complicated functions of the air-off system characteristics. 

There are several notable examples in aeroelastic literature where this category of 
behavior has been produced by analysis or noted in experiment. More notably, however, 
there is a century of experimentation in which this phenomenon has not been observed. 
The current work utilizes a very simple system. In extension to more complicated 
systems, the phenomenon may change or simply be more difficult to observe. This 
category of system behavior has not been widely predicted by analysis and generally the 
dynamic mode behavior is very damped, both factors making it difficult to locate in an 
experimental setting. 

Knowledge of the subcritical, supercritical and noncritical mode behavior is an asset for 
many reasons. Understanding the fundamental physics of a system is a good thing all by 
itself, of course. However, in addition, two practical reasons to have this knowledge 
come immediately to mind. Test techniques for predicting divergence onset have 
included frequency tracking of dynamic modes- divergence onset being indicated by the 
nearness of this frequency to zero. For configurations where divergence occurs without a 
structural dynamic mode losing its oscillatory nature this technique would not alert one to 
the onset of divergence. A second practical reason to understand the noncritical mode 
behavior is related to control applications. As active control of aeroelastic responses 
becomes more commonplace, it becomes more vital to understand the behavior of system 
modes which are noncritical for the uncontrolled or open loop system. Control law 
designs which are model-based rely on modal knowledge of system characteristics, not 
simply stability. 

The body of this paper first presents a discussion of aeroelastic analysis methods and a 
historical perspective of programs which studied the divergence phenomenon. In the 
background material, a discussion is presented of past research that studied divergence 
mechanisms. The present analysis method and results are then presented and discussed. 
These results include a discussion and examples of the parametric database that 
delineates regions where the dynamic mode behavior at divergence changes. Detailed 
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analytical results are presented for one configuration of the wind tunnel model design. 
These results include stability analysis and study of the eigenvectors. Brief results are 
presented for the two additional wind tunnel model configurations. Analytical results are 
also presented for a set of parameters which produces divergence after the structural 
dynamic originated mode has become non-oscillatory. The experiment is next described. 
The model design process is summarized, as well as the hardware employed. The 
experimental techniques and data reduction methods are addressed. The results of the 
experiment are presented for the three configurations tested. The data are presented and 
discussed in the following order: determination of the divergence condition, subcritical 
techniques for predicting divergence onset, system behavior at divergence, and subcritical 
modal characteristics. Analytical and experimental results are then compared in terms of 
the divergence dynamic pressure and modal characteristics. The body of this work 
concludes with a discussion, summary of conclusions and suggestions for future 
directions. 
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CHAPTER ONE 


BACKGROUND 

Aeroelasticity is concerned with problems in which there is substantial interaction among 
the aerodynamic, inertial and structural forces of an object. When a body moves through 
the atmosphere, or when a body is placed in a wind tunnel, aerodynamic forces act over 
its surface. If the body is deformed, there is a change in the magnitude and distribution 
of these surface forces. This redistribution causes additional deformations; the result is 
an interactive feedback loop between aerodynamic loads and aircraft deflections. 1 Static 
aeroelastic behavior is generally considered to be a study of the mutual interaction 
between static aerodynamics and the stiffness, but not the inertia, of an elastic structure. 

Background material is presented on several topics which unite in this work. Methods 
which have been used to examine aeroelastic stability are discussed first. A historical 
look at programs which have studied aeroelastic divergence is then presented. Focusing 
on divergence by virtue of an aerodynamic-originated root, as distinct from a root of 
structural origin, then follows. 


Aeroelastic Stability Methods 

Methods for analyzing the stability of an aeroelastic system set the foundation for the 
work to be presented. The classical methods of solving for stability of the aeroelastic 
equations are the p-method, the k-method and the p-k method. Each of these methods, 
which go by several names, will be discussed in the following paragraphs. In the 
following methods, the non-dimensional Laplace operator, or differential operator, is 
denoted p. In addition to these well-established methods, the g-method will also be 
briefly discussed. 

The linearized equations of motion for a flexible aircraft contain unsteady aerodynamic 
terms, which depend on the Mach number, M, and the reduced frequency, k. For all but 
the simplest aerodynamic theories, the exact aerodynamic coefficients which are 
dependent on M and k, have not been developed in the form of algebraic functions. As a 
result, aerodynamic coefficients are often computed for each desired Mach number for a 
set of predetermined values of reduced frequency. 


1 Weisshaar, Terrence B., Fundamentals of Static Aeroelasticity; Dowell. Earl H., Edward F. Crawley, 
Howard C. Curtiss Jr, David A. Peters, Robert H. Scanlan and Fernando Sisto. A Modern Course in 
Aeroelasticity; Raymond L. Bisplinghoff and Holt Ashley, Principles of Aeroelasticity 
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The k-method is also known as the V-g method or the American method of flutter 
solution to determine the aeroelastic stability ot a system. Many aerodynamic 
formulations, such as the Doublet Lattice method, lead to aerodynamic matrices which 
are only valid for harmonic motion, p=ik. Using these simple harmonic loads, and 
introducing an artificial structural damping factor, complex roots are obtained from the 
equations. There are several well-known disadvantages to the k-method. The complex 
eigenvalues obtained do not represent the actual damping or frequency of the system 
modes except for neutrally stable roots, where the damping is zero. Many solutions are 
required to obtain “matched-point” flutter boundaries. For a given airspeed, several 
solutions with different frequencies may occur. Information regarding non-critical 
conditions and eigenvalues is only qualitative. 

The p-k method, sometimes referred to as the “British Method” or as Hassig’s modified 
version of the Frazer and Duncan method, attempts to improve upon the k-method by 
allowing the reduced frequency to be> complex. In 1971, in discussing the p-k method, 
Hassig wrote, “It is generally conceded that it is desirable to formulate and solve the 
flutter equation such that the solution leads to a vqlue for the rate of decay. Ideally, this 
requires the formulation of the unsteady aerodynamics matrix as a function of the 
complex variable p. When one wants to work with exact theoretical aerodynamics one 
must work with a formulation for harmonic motion and devise approximate methods to 
determine the rate of decay.” The equations of motion are written in a form indicating 
that the aerodynamic matrix is available only for harmonic motion. The eigenvalues of 
this approximate system can be solved, producing complex roots. The aerodynamics are 
then recomputed using the frequency that resulted from the eigenvalue computation. The 
equations of motion in the p-k method are solved in an iterative fashion so that the 
assumed value of k converges to the computed value of the imaginary part of a pre- 
selected eigenvalue. The iterations are repeated, for a single mode at a time, until all the 
modes have achieved convergence. There are several disadvantages to the p-k method. 
While the results for the flutter condition are shown to be quite good, the eigenvalues of 
damped modes are only approximate. The calculated damping is only good for low 
levels of damping. Another disadvantage of the p-k method is the requirement to track 
the eigenvalues of the system as the velocity or dynamic pressure is increased. For an n- 
degree of freedom system, as each mode is tracked, the equations produce n eigenvalues. 
Selection of the proper root is vital to the success of the method. 

The p-method is the simplest method to understand, but perhaps the most difficult to 
apply. Utilizing the p-method means simply solving for the complex eigenvalues of the 
governing equations. Bisplinghoff and Ashley 2 comment on the process of finding 
eigenvalues of an aeroelastic system to determine stability: “The system consisting of a 
typical section in an airstream possesses dynamic eigenvalues. The critical (instability) 
condition is defined to occur at the lowest speed ... at which the damping ratio of any 
aeroelastic mode passes through zero. Mathematically they consist of values of the 
complex (non-dimensional Laplace) variable, p, which cause the determinant ot the 


2 Raymond L. Bisplinghoff and Holt Ashley, Principles of Aeroelasticity 
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(nondimensionalized coupled aeroelastic system) matrix to vanish. Because (the 
aerodynamic velocity potentials) are transcendental functions (of p), there exists in theory 
an infinity of such roots. So, with the air included, (the aeroelastic system) has an infinite 
number of degrees of freedom. . . . The most logical way of studying the dynamic 
aeroelastic stability of a structure. . . could seem to be to calculate the root locus of p as a 
function of airspeed and altitude. In engineering practice, however, this has not been the 
customary approach, as ... more data are available on airloads resulting from simple 
harmonic motion.” The p-method avoids the iteration process by using explicit 
expressions for the aerodynamics. If the aerodynamics can be expressed as a sufficiently 
simple function of p, the aeroelastic equations define a polynomial in p. The main 
difficulty with the p-method lies in the derivation of appropriate aerodynamic 
expressions. The p-method has been used with quasi-steady aerodynamics represented 
by a first order differential equation, which ignores any effect of the wake. Two methods 
of approximating higher order aerodynamic theories are the Pade method by Vepa 3 and 
the minimum state method by Karpel and Hoadley 4 . These methods apply rational 
polynomial fits to the tabular values of the complex aerodynamic coefficients which were 
derived for oscillatory motion. Another method of approximation was developed by 
Nissim 5 , in which a second order complex coefficient fit is used rather than a rational 
function approximation. 

A damping perturbation method, named the g-method, has recently been developed by 
Chen 6 . This is a generalization of the k-method and the p-k method. The basic 
assumption is that a first order Taylor series approximation in g can be developed. Chen 
utilizes analytical continuation to replace the derivative with respect to g with a 
derivative with respect to reduced frequency. He states that this substitution is valid in 
the complete p-domain except along the negative real axis in subsonic flow. The g- 
method produces results that agree well with well-established methods. His method also 
yields some results in which aerodynamic lag divergence is illustrated. He explicitly 
points out that damping results of the p-k method are valid for where damping or reduced 
frequency are zero or where the change (derivative) in the aerodynamics with respect to 
reduced frequency is zero. 


Divergence, down through the ages 

New flight vehicle concepts often invigorate the study of aeroelasticity, as new types of 
interactions are anticipated or observed. Examining the literature on divergence, several 


; Vepa, R., On the use of Pade Approximants to Represent Unsteady Aerodynamic Loads for Arbitrarily 
Small Motions of Wings 

4 Karpel, Mordechay, and Sherwood Tiffany Hoadley, Physically Weighted Approximations of Unsteady 
Aerodynamic Forces Using the Minimum-State Method 

5 Nissim, E., Flutter Analysis Using a New Complex p-Method 

6 Chen, P.C. A Damping Perturbation Method for Flutter Solution: The g-Method. 
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programs or research areas stand out as having inspired the study of aeroelastic 
divergence. 

Consideration of various plantorms sparked work in static divergence at the National 
Advisory Committee for Aeronautics (NACA) in the 1940’s and 1950 s. Diederich and 
Budiansky 7 conducted analytical and experimental research into swept and tapered wings 
which resulted in charts and approximate formulae for estimating wing divergence tor 
various configurations. This work demonstrated among other things, the dramatic 
decrease in the divergence speed for a wing that was swept forward. The analytical work 
investigating different planforms was continued by Diederich and Foss , as they produced 
an analytical method of calculating divergence of wings with various planforms including 
delta wings. A summary of the NACA analytical efforts is provided by Diederich 1 . 


The 1970’s and 1980’s brought a resurgence in the study of static aeroelastic effects. 
Advances in composite materials lead to a reconsideration of the forward-swept wing 
concept, which had been previously dismissed due to divergence. A study by Krone 
showed clearly that “the detrimental effect of divergence on forward swept airfoils can be 
successfully controlled. By tailoring the composite layer thickness distributions and 
orientations a design can be obtained that produces optimum stiffness and strength 
characteristics ... with little fear of suffering the weight penalties that have previously 
been caused by the divergence phenomenon." Static aeroelastic characteristics of 
forward swept wings were investigated by many aeroelasticians. To note a few, 
Weisshaar 1 1 discussed forward-swept wing divergence from a fundamental concepts 
point-of-view, and Blair 12 performed wind tunnel experiments which demonstrated the 
fundamental relationships among sing sweep, composite fiber orientation and divergence 
speed. An experimental study of the static aeroelastic divergence of forward-swept 
wings was conducted in the NASA Langley Transonic Dynamic Tunnel by Ricketts and 
Doggett 13 . Flat plate models with varying geometry were tested. Six subcritical response 
testing techniques were formulated and evaluated at transonic speeds for accuracy in 
predicting static divergence. Ricketts and Doggett concluded that, “in general, the static 
methods seemed to consistently give better quality data than the dynamic methods.” As 
pointed out by Doggett and Ricketts, dynamic methods of divergence prediction produce 
inferior results to those produced by static methods. The current work addresses the issue 
of why this is true. The dynamic behavior is governed not only by the stiffness (static) 
properties but by the inertial properties. Dynamic methods work only if a complex mode 


7 Diederich. Franklin W., and Bernhard Budiansky. Divergence of Swept Wings 

s Diederich, Franklin W., and Kenneth A. Foss. Static Aeroelastic Phenomena ofM-, W-. and -Wings. 

11 Diederich. Franklin W., Divergence of Delta and Swept Surfaces in the Transonic and Supersonic Speed 
Ranges 

10 Krone, Norris J., Jr., Divergence Elimination with Advanced Composites . 

11 Weisshaar, Terrence B. Forward Swept Wing Static Aeroelasticity. 

12 Blair. Maxwell, Wind Tunnel Experiments on the Divergence of Swept Wings with Composite Structures 

13 Rodney H. Ricketts, and Robert V. Doggett, Jr, Wind-tunnel Experiments on Divergence of Forward - 
Swept Wings . 
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of the system becomes real. A survey article by Shirk, Hertz and Weisshaar 14 provides 
an extensive reference list for other work on this subject. 

Two forward swept wing airplanes were designed, built and flight tested. “Grumman 
Aircraft Corporation built two X-29's. Phase 1 of the project, using aircraft No. 1, was 
flown from December 1984 to 1988 and investigated handling qualities, performance, 
and systems integration. Phase 2 of the X-29 program involved aircraft No. 2 and studied 
the high angle of attack characteristics and military utility of the X-29. 15 “ In the 
development of the flight vehicle concept, dynamic analysis and wind tunnel testing of a 
free-free configuration was performed by Miller, Wykes and Brosnan 16 . Their analyses 
revealed a different type of instability. The phenomenon involved a coupling between 
the wing divergence mode and the aircraft short period mode, termed rigid body/ wing 
bending flutter. The analytical results showed that the wing response was completely 
different from the cantilevered case. While the divergence of the forward swept wing 
flight vehicle was controlled by aeroelastic tailoring, the coupling of the wing 
divergence-prone mode and the rigid body motion was controlled by enhancement of the 
stability augmentation system (SAS). One of the research objectives for the X-29 flight 
test program became correlating flight data with the predicted structural stability and 
determination of the aeroservoelastic stability margins. 17 In testing the X-29, methods of 
divergence prediction as applied to flight tests were investigated. Schuster and Lokos 
focused on applying the Southwell method to flight test data. Consideration of potential 
errors lead to a more conservative pace in envelope expansion than might otherwise have 
been required. 

The 1990’s topic in divergence centered around the National Aerospace Plane (NASP). 
Researchers working on this program examined divergence of all-moveable surfaces, 
which were representative of wing configurations under consideration. Experimental 
data was obtained in a supersonic test conducted in the Unitary Plan Wind Tunnel at the 
NASA Langley Research Center 18 . The wing models had low aspect ratios and highly 
swept leading edges. The wings were attached by a single-pivot mechanism along the 
wing root. The supersonic divergence was predicted to be primarily dependent on the 
first wing pitch mode. Two subcritical response instability prediction techniques were 
used: the static Southwell method and the dynamic frequency tracking method. The 
improved Southwell method uses the change in slope of load-versus-angle of attack 
measurements as dynamic pressure is increased to predict divergence conditions. 

Accurate predictions were not obtained for this wind tunnel model using measurements 
from the strain gauge bridges on the pitch stiffness elements. The frequency of the wing 


14 Shirk, M.H., T.J. Hertz, and T.B. Weisshaar. Aeroelastic Tailoring- Theory, Practice, and Promise. 

15 http://www.dfrc.nasa.gov/gallery/photo/X-29 

16 Miller, Gerald D., John H. Wykes and Michael J. Brosnan. Rigid Body-Structural Mode Coupling on a 
Forward Swept Wing Aircraft 

17 Sefic, Walter J., and Cleo M. Maxwell, X-29 A Technology Demonstrator Flight Test Program Overview 

18 Stanley R. Cole, James R. Florance, Lee B. Thompson, Charles V. Spain and Ellen P. Bullock, 
Supersonic Aeroelastic Instability Results fora N ASP-like Wing Model. 
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pitch mode was tracked as the dynamic pressure was incrementally increased in the wind 
tunnel, and the model vibration frequency as it approached zero at the divergence 
condition was monitored. Extrapolation of subcritcal measured frequencies was then 
performed. This method was successfully used during the wind tunnel test to extrapolate 
divergence dynamic pressure and as guidance in anticipating actual divergence 
instabilities. 

It is interesting to note that the experimental data show that divergence occurred at 
dynamic pressures only 3-6 percent beyond the dynamic pressure at which the frequency 
reaches a value that is fifty percent of the wind off natural frequency. As described by 
the authors, there is a large change in frequency that is determined only by observing the 
dynamic pressure of the instability itself. For two of the Mach numbers for the nominal 
stiffness configuration, the nearest subcritical data indicates that the mode has a 
frequency that is at 50% of its wind off value. For the other two cases, however, there is 
subcritical data which show the pitch mode frequency has dropped much further just 
prior to divergence. 


N on-traditional divergence 

Static divergence that occurs without a structural dynamic mode losing its oscillatory 
nature and becoming static is central to the current work. A very interesting body of 
work on this subject exists in the literature. Studies will be discussed which were 
performed by Rodden and various co-authors, Edwards, Dashcund and Martin and 
Watkins. Already mentioned in the discussion of stability methods is the work of Chen, 
in which aerodynamic lag divergence was also found. 

From 1969 to 1994, publications by Rodden and various coauthors present analytical 
results which demonstrate aerodynamic lag divergence and provide a method for 
calculating the true damping of non-critical modes. In the first of these articles, Rodden 
and Stahl 19 performed aeroelastic stability analysis utilizing the p-method. A transient 
formulation of the flutter and divergence problems was presented using aerodynamic 
strip theory and an exponential approximation (the W.P. Jones approximation) to the 
Wagner function. The limitation of the method presented in the referenced work is in the 
aerodynamic strip theory approximation. 

A cantilevered wing with 5 structural modes was analyzed. The divergence velocity was 
found and agreed very well with static calculations. Using the p-method, they 
determined the subcritical frequencies and dampings. They discovered that tracking the 
frequency of the mechanical modes did not produce the instability. Tracking the 
aerodynamic lag roots, which are explicitly present due to the W.P. Jones approximation, 
produced the divergence instability. Their results show the first mode frequency curve 
decreased rapidly at speeds slightly higher than the divergence speed. The frequency 


19 William P. Rodden and Bernhard Slahl. A Strip Method for Prediction of Damping, in Subsonic Wind 
Tunnel and Flight Flutter Tests. 
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went to zero in a small range of velocity and then increased rapidly. Their results are also 
show that the subcritical damping values differ significantly from the artificial damping 
predicted by another stability analysis method, except for the instances where the modes 
behave in a nearly simple harmonic fashion. 

Additional work published by Rodden, Harder and Bellinger' 0 compared the p-method 
results from the above work to results utilizing a p-k solution. While the divergence 
velocity predicted by the two methods agreed, this publication indicated that the 
divergence mechanisms predicted disagreed. The p-method indicated that the mechanism 
was aerodynamic lag divergence; the p-k method predicted divergence of the first 
bending mode. The matter was revisited' and the p-k results reinterpreted. The 
transition from the bending mode to the aerodynamic lag root was then recognized, 
bringing the predicted divergence mechanisms into agreement- an aerodynamic lag 
divergence. 

The previously cited references analyzed a cantilevered wing. This work was extended to 
include a vehicle plunge mode". Applying the aerodynamic approximations and p- 
method analysis technique predicted instability of an aerodynamic lag root. In the 
unrestrained system the instability is oscillatory as the unstable aerodynamic mode is 
coupled with the vehicle plunge freedom. 

The cantilevered configuration was revisited by Rodden and Johnson in 1994' . The 
subsonic Doublet-Lattice aerodynamic method was employed in a p-k solution 
procedure. This analysis shows no aerodynamic lag root present in the divergence 
mechanism. They commented that the first bending frequency moved smoothly to zero 
frequency. They assert in this publication that the demonstrated discontinuous behavior 
of the eigenvalues is not due to a physical phenomenon, but due to the change in the 
definition of damping when a root becomes real. The authors do not comment in this 
publication on anticipated inaccuracies due to aerodynamic modeling or inexactness of 
the p-k solution for predicting subcritical characteristics. 

Aerodynamic mode divergence was also illustrated in an analytical study by Edwards . 
In this work, he discussed aerodynamic modeling in depth and notes some of the 
shortcomings of different aeroelastic stability methods. 

Edwards began with a solution of the linearized potential equation for the case of two- 
dimensional airfoils undergoing simple harmonic motion in incompressible flow, 
published by Theodorsen. He next extended the derivation to arbitrary motion or 

20 William P. Rodden, R.L. Harder and E. Dean Bellinger, Aeroelastic Addition to NASTRAN. 

21 William P. Rodden and E. Dean Bellinger, Aerodynamic Lag Functions , Divergence and the British 

Flutter Method. . 

22 William P. Rodden and E. Dean Bellinger, Unrestrained Aeroelastic divergence in a Dynamic Stability 

Analysis. 

22 William P. Rodden and Erwin H. Johnson. MSC/NASTRAN Aeroelastic Analysis User s Guide. 

24 John E. Edwards, Unsteady Aerodynamic Modeling and Active Aeroelastic Control. 
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complex values of reduced frequency. He called this the generalized Theodorsen 
aerodynamic representation. The expressions for the aerodynamic loads, lift and pitching 
moment, were incorporated into typical section equations of motion. The stability of the 
aeroelastic system was investigated using the p-method. The advantage of the p-method 
is that the exact roots and the degree of stability of the system are determined, to the 
extent of the accuracy of the aerodynamic representation. Edwards’ presented a 
derivation and analytical results which produce an aerodynamic mode at divergence. 

This mode produces the motion of the diverging airfoil and occurs in addition to the 
structural poles. Inherent in his work is the realization that the aerodynamic equations 
are not constant coefficient equations. The fundamental theorem of algebra 2 * states that 
an nth order polynomial equation is guaranteed to have exactly n roots for the case of 
constant coefficients. In the case of the Theodorsen aerodynamic representation, the 
coefficients of the governing polynomial are not constants and as such, no guarantee as to 
the number of roots can be asserted. 

The occurrence of this divergence mode was studied by locating the poles of the system 
in the complex plane. Both the exact system model and a Pade approximate model were 
used to locate these poles. The diverence speed was indicated for the exact model by the 
emergence of an additional real pole on the positive real axis. The Pade model contained 
an eigenvalue which migrated from the stable negative real axis into the unstable positive 
real axis. Both results produced the same value for divergence speed. 


Divergence in the case of a wing instead of a typical section was investigated by 
Dashcund* 6 , the distinction being that his model had wing modes, not rigid pitch and 
plunge degrees of freedom. In the course of performing a flutter suppression study using 
active control, Dashcund discovered that for his configuration a non-structural originated 
root diverged. Divergence is mainly addressed by the analytical portion of the work. 

The equations of motion were generated employing a Rayleigh-Ritz energy method. The 
included modes were beam bending modes and rod twisting modes. The aerodynamics 
were lull unsteady 2-D strip theory. Dashcund writes: “Modeling the unsteady aero in 
the Laplace domain in terms of an irrational, exact representation of the generalized 
Theodorson function shows the presence of additional stability roots which are not 
associated with the structural modes of the system nor with the feedback compensation or 
control surface actuator dynamics. The existence of these additional aerodynamic system 
roots, which includes the divergence root, is confirmed by the qualitatively good 
agreement between predicted and experimental divergence boundaries for the active 
flutter controlled wing.” These results provide an indication that aerodynamic-based 
divergence can exist on wings. A recommendation presented in this work is that a 
frequency tracking of all modes needs to be performed experimentally. This would show 
that the still-existing structural dynamic modes are stable, while simultaneously 
observing the divergent mode. 


* Johnson, R.E, and Fred L. Kiokemeisler. Calculus 



Martin and Watkins 26 present and discuss analytical and experimental data for delta 
wings. The data presented is of a summary nature. Only the divergence conditions are 
given; no frequency information is provided. In their discussion of the transonic test 
data, however, they say, “The divergence dynamic pressures were very sharply defined 
and were marked by one or two large excursions of the tip of the model just prior to 
divergence. ... The model motion, when divergence was reached, was quite rapid and 
the deflection quickly increased until the model was bent beyond 90 degrees to the 
airflow.” Their comments lead to speculation regarding the nature of the divergence 
mechanism that they observed. 


2(1 D. J. Martin and C. E. Watkins. Transonic and Supersonic Divergence Characteristics of Low-aspect- 
ratio Wings and Controls. 


12 


CHAPTER TWO 


ANALYSIS 

A standard procedure for solving a structural dynamic problem is to employ 
eigenanalysis to calculate the structural dynamic eigenvalues and eigenmodes. Recently, 
this eigenvalue/eigenmode procedure has been extended to unsteady aerodynamics, and 
to coupled aeroelastic equations'. 

In computational fluid dynamics, CFD, there are two approximations that are typically 
employed. One is the construction of a computational grid, which determines the limits 
of spatial resolution of the computational model. The second is the approximation of an 
infinite fluid domain by a finite domain. It is a principal purpose of the present 
discussion to demonstrate that the computational grid not only determines the spatial 
resolution obtainable by the CFD model, but also the frequency or temporal resolution 
that can be obtained. Also, as will be shown, the finiteness of the computational domain 
determines the resolution of the eigenvalue distribution for a CFD model. Both of these 
observations have important ramifications for assessing the CFD model and its ability to 
provide an adequate approximation to the original fluid model on which it is founded. To 
these ends, a finite-wake, time-domain, discretized vortex lattice aerodynamic model has 
been utilized. 

Results of aerodynamic parametric variations are presented, as well as detailed discussion 
of the trends produced by these systematic variations. The discussion includes the 
parametric effects on both the discrete- and continuous-time aerodynamic eigenvalues. 
These studies give insights into aerodynamic modeling in the discrete time domain 
including how one may construct reduced order aerodynamic models using the dominant 
aerodynamic modes. 

The aerodynamic model was also combined with time-domain discretized structural 
dynamic equations to examine the aeroelastic behavior of a typical section. Aeroelastic 
response is also discussed in terms of eigenanalysis results. Aeroelastic stability analyses 
generally focus on the migration of the eigenvalues as a function of the velocity or other 
flow parameter. Indeed, much flutter analysis in practice today uses at best only an 
approximation to the true aeroelastic eigenvalues. Here, the true eigenvalues are found 
for all aeroelastic modes without iteration. This enables an examination of the subcritical 
modal characteristics of the system as well as the behavior the noncritical modes at 
instability. 


1 Dowell, Earl H.. Kenneth C. Hall, and Michael C. Romanowski, Eigenmode Analysis in Unsteady 
Aerodynamics: Reduced Order Models: Hall. Kenneth C.. Eigenanalysis of Unsteady Flows about 
Airfoils, Cascades and Wings. 
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The eigenvalue migration is observed and the eigenvectors are investigated to determine 
the relative participation of the aerodynamics and the structural dynamics. 
Nondimensional parametric variations were performed to investigate the changing 
character in the dynamic modes as the system diverges. Using this database, 
configurations were identified and analyzed which exhibited different types of dynamic 
mode behavior as the typical section became statically unstable. Specific emphasis is 
paid to configurations where the dynamic mode persists at a nonzero frequency as the 
system destabilizes. 

The eigenvectors associated with the dynamic and divergent modes are also studied. 
They provide a wealth of information and can supplant or supplement the eigenvalues in 
providing stability information. They are studied from the standpoint of their own modal 
content as well as their phase relationship to other eigenvectors of the system. They are 
utilized to identify the relative importance of the structure and the aerodynamics in a 
given aeroelastic mode’s behavior. 

Aerodynamic Studies 
Aerodynamic Modeling 


A Vortex Lattice solution to Laplace’s equation for incompressible two-dimensional flow 
is utilized in this study. The flow over an airfoil with a certain number of vortex 
elements on the airfoil and in the wake is now considered. The airfoil is modeled as a 
two-dimensional flat plate. The airfoil and the wake are divided into segments, referred 
to as aerodynamic elements. Vortex lattice aerodynamics are generated by placing 
vortices of strengths to be determined at points on the airfoil and in the wake. 

Collocation or control points, usually located aft of the vortex locations, are points where 
the boundary conditions must be satisfied. Typical placement is for the vortices to be 
located at the 14-chord points of the aerodynamic elements. The collocation points are 
typically placed at the 3 4-chord locations of the elements. 

The governing equations are presented by Hall" and shown in detail in Appendix A; they 
are briefly summarized here. There are 3 basic relationships, described in the following 
paragraph, which are combined to form a matrix equation for the vortex strength. 
Equation 1, where n and n+1 denote the next and the current discrete time sample. T is a 
vector of vorticities and w is a vector of downwashes at each of the collocation points. 
The number of elements on the wing is denoted M, while the total number of elements is 
denoted N. 

[Afrf +] + [Bfrf = Vf +I Equation I 


2 Hall. Kenneth C„ Eigeuanalysis of Unsteady Flows about Airfoils. Cascades and Wings. 
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Three basic relationships determine the contents of the A and B matrices seen in Equation 
I . These represent N equations with N variables. The first of the three basic 
relationships equates the velocity induced by the discrete vortices at the collocation 
points to the downwash caused by the airfoil’s motion. This relationship accounts for M 
rows within the matrices where M is the number of spatial grid points on the airfoil or 
wing. Applying Kelvin’s theorem generates a second basic relationship utilized in 
deriving the matrix equations. Unsteady vorticity is shed into the wake; its strength is 
proportional to the time rate of change of circulation about the airfoil. The time step is 
taken to be equal to the time it takes the vorticity to convect from one vortex station to 
the next. This relationship accounts for the (M+l ) row of the matrix equations. Once the 
vorticity has been shed into the wake, it is convected in the wake at the freestream 
velocity. This is the third basic relationship which appears in Equation 1 as rows (M+2) 
through (N-l). Vorticity convection also provides the final, Nth, row of the matrix 
equations. Because the wake is modeled with a finite length, the last vortex element must 
be treated specially. “Otherwise, the starting vortex would disappear abruptly when it 
reached the end of the computational wake, producing a discontinuous change in the 
induced wash at the airfoil. To alleviate this difficulty, ... the vorticity is allowed to 
dissipate smoothly by using a relaxation factor.” * 

The formulation and analysis of the aerodynamic model progresses in the following 
manner. Discrete, time-marching equations are written as shown in Equation 1. Once 
these equations are written, they inherently contain the approximations of the finite w ake 
and the discretization. A discrete Fourier transformation is performed on the unforced 
equations, producing the z-plane representation, Equation 2. 

ZT() = (- A - ' Equation 2 


The discrete time eigenvalues, z, and the eigenvectors, r„, are extracted from these 
equations. These provide insight into the behavior of the aerodynamic model and also 
provide a method for constructing a reduced order model. These eigenvalues are then 
converted to the continuous time domain, X-plane, through a zero order hold 
transformation, Equation 3. 

2 = jgg(g) 


At 


Equation 3 



Baseline Configuration 


As the first of several numerical examples, the flow over an airfoil with 20 vortex 
elements on the airfoil and 180 elements in the wake, equally spaced, is now considered. 
This will be referred to as the baseline case. The (finite) length of the wake thus extends 
9 chord lengths. The eigenvalues and eigenmodes of the flow can be computed by 
established methods. Because there are 200 elements in the model, 200 eigenvalues 
result. 

The discrete time (z-plane) eigenvalues, extracted from Equation 2, approximately form a 
circle centered at the origin, as shown in Figure 1. In addition to these eigenvalues, there 
are a finite number of eigenvalues at the origin. The number of eigenvalues at the origin 
is equal to the number of segments or grid points on the wing. This conclusion follows 
from examining the rank of the system matrices in equation 1, from the numerical results 
obtained here, and appears to be supported by the results presented in Hall', though it was 
not noted in this previous work. Eigenvalues at the origin in the discrete time domain 
transform to -« in the continuous time domain. 



Real Part 


Figure 1 Eigenvalues for baseline case; discrete time eigenvalues, z 
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The continuous time eigenvalue distribution for the baseline case is shown in Figure 2. 
The real part of the eigenvalue is indicative of the damping and the imaginary part is the 
damped frequency of each fluid eigenmode. Examining the eigenvalues of the 
aerodynamic matrix in the continuous domain produces several observations. The 
continuous domain eigenvalues are discretely spaced and are arranged in “arms" that 
emanate from the origin and reach up and down in the left half plane. Additionally, the 
real parts of the arms asymptotically approach a limiting value. 



Real Part 

Figure 2 Eigenvalues for baseline case; continuous time eigenvalues, X 

The presence of positive aerodynamic damping is evidenced by the arms lying in the left 
half plane. The primary contribution to the damping appears to lie with the overall flow 
field, however, there is additional damping due to the presence of a vorticity relaxation 
factor at the last wake element. The relaxation factor used in the vortex lattice model 
provides energy dissipation in the wake; as the relaxation factor is decreased, more 
energy is dissipated and the aerodynamic damping increases. If the number of 
aerodynamic boxes within the wake is increased, the last box will be a smaller percentage 
of the total wake length and thus, the influence of the relaxation factor will be 
diminished. 
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Parametric Variations 


Three aspects of the aerodynamic modeling significantly impact the eigenvalue 
distribution: the size of the aerodynamic elements, the number of these elements that lie 
in the wake, and the length of the wake. The three aerodynamic configurations, detailed 
in Table 1, compared against each other two at a time, produce the three comparison 
cases, which are organized in Table 2 and discussed next. 


Aero Config 
No. 

Airfoil 

Wake i 

No. of 
elements 

Normalized 
element size 


No. of 
elements 

Normalized 
element size 

Normlized 
wake length 

1 

(Baseline) 

20 

1 

1 

180 

1 

1 

2 

20 

1 

1 

360 

1 

2 

3 

40 

Vi 

1 

360 

Vi 

1 


Table 1 Aerodynamic Configurations 


Comparison 
Case No. 

Aerodynamie Configurations 
Compared 

Parametric Variation 

Quantity Held Constant 


1 

2 

3 



I 


X 

X 

Size of aerodynamic 
elements in wake 

Number of aerodynamic 
elements in wake 

II 

X 

X 


Number of aerodynamic 
elements in wake 

Size of aerodynamic 
elements in wake 

III 

X 


X 

Size and number of 
elements in wake 

Length of wake 


Table 2 Comparison Cases for Parametric Variations 


The three comparison cases are discussed in terms of their discrete time eigenvalue 
distributions (z- values), their discrete-to-continuous time domain transformations (z- 
transformations) and their continuous time eigenvalue distributions (A,-values). 

Comparison case I compares aerodynamic configurations 2 and 3, examining the effects 
of varying the size of the aerodynamic elements while maintaining the number of 
elements which lie in the wake. Because the number of wake elements remains fixed, 
configuration #2 has a wake that is twice the length of the wake in configuration #3 and 
elements which are twice as large. Although not shown, the discrete time eigenvalue 
patterns for configurations 2 and 3 are identical because the number of elements in each 
wake is identical. However, changing the size of the aerodynamic elements changes the 
transformation, which must be applied to convert the discrete time system to continuous 
time. This difference in transformation produces the change in continuous domain 
eigenvalues, as illustrated in Figure 3. 
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Figure 3 Case I: Influence of varying the size of the aerodynamic elements. Continuous time 
eigenvalues, X 


It is easily shown that the frequency of each eigenvalue scales linearly with the 
aerodynamic element size. The maximum frequency of the arms can be determined a 
priori by utilizing Shannon’s sampling theorem. The aerodynamic eigenfrequencies are 
bounded from discrete time considerations similar to those that predetermine the discrete 
Fourier transform frequencies 3 . The maximum frequency, ft), that can be resolved would 
have 1 cycle spanning two aerodynamic panels. Using the velocity to relate the spatial 
and temporal sample sizes. Equation 4, leads to maximum frequency that can be resolved. 
Equation 5. 

Ax 

\J = Equation 4 

At 


max(<w)= 


nU 

Ax 


Equation 5 


* 3 Hardin, J.C., Introduction to Time Series Analysis. 
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Thus, changing the aerodynamic element size changes the frequencies of the 
aerodynamic eigenvalues. As the size of the elements becomes infinitesimal, it is 
speculated that the eigenvalue arms will cover the frequency range from ±» . 

It should be noted in studying Case I that the number of eigenvalues has remained 
constant in going from configuration 2 to configuration 3, while the frequency range has 
doubled. Thus, the density of the eigenvalues has halved. The implications of this will 
be further discussed in studying Case III. 

Comparison case II compares aerodynamic configurations 1 and 2 and examines the 
effect of varying the number of aerodynamic elements in the wake while holding their 
size constant. The number of aerodynamic elements in the wake determines the number 
of discrete time eigenvalues comprising the pseudo-circular pattern. As more elements 
are placed in the wake, the more crowded pattern expands outward towards the unit circle 
and the damping of each aerodynamic mode. As the element size decreases, the radius of 
the pseudo-circular pattern asymptotically approaches 1. In discrete time eigenvalue 
analysis, an eigenvalue lying on the unit circle represents a neutrally stable system. In 
the continuous time domain, the imaginary axis is the line of demarcation for stability. It 
is thus anticipated that the additional boxes in the wake force the “arms” of the 
continuous time eigenvalues closer to the imaginary axis. Figure 4 bears this out. As 
more elements are added to the wake, the closer the aerodynamic roots get to those 
associated with simple harmonic motion. Thus, changing the number of aerodynamic 
elements in the wake changes the damping of the aerodynamic eigenvalues. As the 
number of elements goes to infinity, it is speculated that the arms will move to the 
imaginary axis. 
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Figure 4 Case II: Influence of varying the number of aerodynamic elements in the wake. 

Continuous time eigenvalues, X 

It should be noted in studying Case II, as the wake length is increased, leaving the size of 
the aerodynamic elements constant, the frequency range of the continuous time 
eigenvalues remains constant. The number of aerodynamic elements determines the 
maximum frequency. Doubling the number of elements in the wake means doubling the 
number of eigenvalues on the “arms.” Twice as many eigenvalues reside in arms of the 
same length. Hence, the continuous time eigenvalue distribution has become denser. 

Comparison case III compares aerodynamic configurations 1 and 3 and examines the 
effects of varying simultaneously and in inverse proportion, the number and length of 
aerodynamic elements in the wake, such that the wake length remains constant. The 
expected trends for the behavior of the arms of the continuous time eigenvalues are 
difficult to predict because, in going from configuration 1 to configuration 3 there are 
multiple tendencies: increasing the number of elements tends to move the arms closer to 
the imaginary axis; decreasing element size tends to extend the frequency range of the 
arms. The combined result on the continuous time eigenvalues, shown in Figure 5, is that 
the arms of the eigenvalues lie approximately the same distance from the imaginary axis, 
while the frequency range of configuration 3 is twice that of configuration 1. This 
corresponds to the effects of smaller element size of configuration 3. Thus, the spacing 
of the eigenvalues is approximately constant between the two analysis runs. 
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Figure 5 Case III: Influence of simultaneously varying the size and number of aerodynamic elements 
in the wake, maintaining a constant wake length. Continuous time eigenvalues, X 


An approximate formula for eigenvalue spacing is derived using the frequency range and 
the number of eigenvalues. The maximum frequency was found using Equation 5. 
Accounting for positive and negative values, the frequency range is twice this. Dividing 
this range by the number of elements or eigenvalues in the wake, and recognizing that the 
element size times the number of elements in the wake is the wake length produces the 
relationship given in Equation 6. 


Ao) = 


2nU 

I-wake 


Equation 6 


The reader may recognize that this is similar to determination of the discrete Fourier 
transformation frequencies, as determined by the length of the time record. The 
eigenvalue spacing is approximate due to the eigenvalues not lying on the imaginary axis, 
that is, due to the discretization-induced damping. For the case of the element size 
becoming infinitesimally small, the formula is exact. 
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Thus, the effect of the finite wake is to produce discretely spaced eigenvalues, instead of 
a continuous line. As the wake length becomes infinite, it is speculated that the arms of 
discretely spaced eigenvalues form continuous lines emanating from the origin. 


Discussion of Aerodynamic Studies 

The study of aerodynamic eigenvalues using the vortex lattice code has led to some basic 
insights. The eigenvalues have been shown to be artifacts of the discretization and the 
finite length wake. 

The effects of discretization are controlled by two independent factors. The size of the 
elements determines the range of frequencies covered by the eigenvalues, while the 
number of elements in the wake drives the damping. Their effects are shown to be 
independent, as one controls the transformation from discrete to continuous time, and the 
other controls the discrete time eigenvalue pattern. The effect of the finite wake is to 
produce discretely spaced eigenvalues, instead of a continuous line. 

The following speculations regarding the limiting cases are offered. As the size of the 
elements becomes infinitesimal, the eigenvalue arms will cover the frequency range from 
±oo. As the number of elements goes to infinity, the arms will move to the imaginary 
axis. As the wake length becomes infinite, the arms of discretely spaced eigenvalues 
form continuous lines emanating from the origin. 

Aerodynamic eigenvalues have been shown to be artifacts of the discretization, which 
exist regardless of the airfoil or wing motion applied to the model. The eigenvalues exist 
even with no airfoil or wing motion. A direct analogy with the feedback control problem 
can be drawn for aeroelastic systems. The poles of the controller exist, even when the 
system is open loop. The system is open loop when the feedback path is cut. Three 
scenarios produce open loop behavior: the sensor information is not provided to the 
control law, the controller output is not applied to the physical system, or the control law 
has a zero gain. The last case is analogous to the aeroelastic feedback scenario when the 
velocity is zero. Just as the poles and zeros of the control law are independent of the 
feedback gain, the poles or eigenvalues of the aerodynamic system are independent of 
velocity. 

It should be noted that this analogy is not be carried further because standard root locus 
rules of migration for increasing gain are not directly applicable to the aeroelastic 
scenario, except with the simplest aerodynamic models. The open loop aerodynamic 
poles are complicated functions of the velocity, which vary with airspeed. 
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Aeroelastic Studies 


The discrete time aerodynamic model can be coupled with a discretized structural 
dynamic model to produce the following time-marching aeroelastic equations of motion, 
which can then be analyzed to determine the behavior of the system. The vector q 
contains the structural dynamic degrees of freedom, the vector f represents the 
aerodynamic loads and the matrices, D, and D 2 describe the coupling between the 
aerodynamic and structural dynamic quantities present in an aeroelastic system. 

D 2 q' l+] + Dfl" + / W+l = 0 Equation 7 

The aerodynamic loads, f, can be expressed in terms of the unsteady vorticities on the 
wing, T. 

y-rt+l = c~,r n+] +C\T n Equation8 

For a system with no external disturbances, the downwash on the airfoil, w, is produced 
by the motion on the airfoil. 


w n =Eq 


n 


Equation 9 


Combining Equation 1, 7, 8, and 9 produces the aeroelastic system equations, Equation 

10 . 


D 2 C 2 
-E A 


n+\ 


+ 


D, C, 

0 B 



0 

0 


Equation 10 


The Typical Section 

The typical section is a structural and aerodynamic idealization where the motion and the 
airflow can be represented as two-dimensional. The airfoil section is considered rigid 
and its permitted motion limited to vertical translation and rotation about a fixed axis. 
Here, the typical section motion has been further limited to permit only rotation. The 
boundary condition or mounting system is such that the structural stiffness is represented 
by a torsional spring. The axis of rotation is termed the elastic axis; its position is 
measured positive aft from the center of pressure. The geometric parameters are 
illustrated in Figure 6. The non-dimensional parameters of interest for a single degree-of- 
freedom typical section are defined in Table 3. 
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Figure 6 Typical section with pitch freedom 


Description 

Parameter Symbol 

Relationship with 
dimensional quantities 

Torsion mode frequency 

C0 U 

"V A* 

Elastic axis location 

e/b 

e/b 

Mass ratio 

ft 

m/(span p air b 2 ) 

Radius of gyration 

r« 

V 

'la/ 

/ 

/ mb “ 

Reduced velocity 

V 

U/(o«b 


Table 3 Nondimensional parameters of aeroelastic system 


Details of the structural dynamic equations are presented in Appendix A. They are 
represented in generic notation in Equation 7. The generalized coordinate vector, q, 
contains only a single degree of freedom, i.e. angle of attack, and its time derivative. 


Stability Analyses 

The stability of the aeroelastic system was analyzed by solving the equations of motion 
for a series of reduced velocities. Eigenanalyses of the discrete time systems were 
performed on each set of equations and the system eigenvalues tracked. The eigenvalues 
were transformed into the continuous time domain using a zero order hold 
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transformation. Stability can be inferred from either the discrete or the continuous time 
root locus. 


A detailed look at the stability analysis for wind tunnel model configuration # 2 is 
presented. The parameters used in this analysis are summarized in Table 4. This 
configuration diverges while the dynamic mode that originated as the structural mode 
persists. 


Description 

Parameter 

Value 

Units 

Semi-chord 

b 

4 

inches 

Span 

span 

21 

inches 

Radius of gyration 

r« 

0.459 


mass ratio 

fi 

51.42 


torsion mode frequency 

Wa 

49.5 

radians/second 

Elastic axis location 

e/b ' 

0.375 


Number of aerodynamic 
elements on the wing 

M 

• 10 


Total aerodynamic 
elements 

N 

100 


Density of air 

Pair 

0 

— 7 — 

si inches/inch 

Aerodynamic relaxation 
factor 

a 

0.996 


Size of aerodynamic 
element 

Ax 

0.8 

inches 


Table 4 Parameter values used in analysis, wind tunnel model configuration #2 


The discrete time root locus is presented in Figure 7. These z-plane plots show the 
imaginary part versus the real part of the eigenvalues. The structural-dynamic-originated 
mode eigenvalue and the aerodynamic-originating eigenvalues, referred to collectively as 
the aeroelastic eigenvalues, migrate as the reduced velocity is increased. Figure 7 
somewhat resembles the plot of the eigenvalues for the uncoupled aerodynamic equations 
which was presented in Figure 1. The complex aerodynamic-originating eigenvalues 
appear relatively undisturbed by the coupling with the structural dynamic equations. In 
addition, the single structural dynamic eigenvalue can be seen near the unit circle, 
indicating that it is more lightly damped than the aerodynamic eigenvalues. It undergoes 
substantial movement with the increase in velocity. 

An instability occurs when an eigenvalue lies outside the unit circle. For this system, this 
is observed on the positive real axis. The axes are expanded to more closely examine the 
behavior near instability, Figure 8. This figure shows the migration of the structural- 
dynamic-originating eigenvalue, and also the interplay with several aerodynamic 
eigenvalues. The lowest complex aerodynamic eigenvalue is clearly influenced, as well 
as the real aerodynamic eigenvalues, one of which becomes unstable. It is difficult to 


26 




further study system behavior from these graphs because each velocity produces 
eigenvalues that essentially belong in different z-planes. This will be discussed in detail 
in a subsequent section of this paper. 



Figure 7 Discrete time root locus for configuration # 2 
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Real Part 

Figure 8 Discrete time root locus for configuration # 2, expanded scale 

The aeroelastic system is converted to the continuous domain by zero order hold 
transformations. The behavior of the continuous time domain eigenvalues is shown in 
Figure 9. For clarity, only the region near the origin is presented. The influence of 
velocity on the aerodynamic eigenvalues is now evident. As in the aerodynamic case 
previously discussed where changing the size of the aerodynamic elements changed the 
transformation from the discrete to continuous time, the same effect is now observed for 
the aeroelastic case as the time step size is changed. Recall that the aerodynamic 
parametric studies were conducted at a fixed velocity. The aerodynamic eigenvalue 
“arms” are stretched with increasing reduced velocity. As velocity increases, the 
individual complex eigenvalues’ frequencies increase at constant damping. 
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Real Part 


Figure 9 Continuous time root locus for configuration #2 

Increasing velocity produces a migration in the structural dynamic mode also. The 
coupled mode that originates as the structural dynamic mode will be referred to here as 
simply the dynamic mode of the system. This mode is a pure structural mode only at 
zero airspeed. For any finite velocity, it and all other modes are strictly speaking 
aeroelastic modes. The lowest reduced velocity for which this system was analyzed was 
0.2. The structural dynamic mode for this nearly-zero velocity is indicated by a solid 
triangle in Figure 9. The root lies at 49.5 rads/second, which agrees with the torsional 
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natural frequency. It is helpful to simultaneously examine Figure 9 and Figure 10 when 
interpreting root migration. Figure 10 shows the real and imaginary parts of each 
eigenvalue plotted versus reduced velocity. 
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Figure 10 Continuous time eigenvalues for configuration # 2 as functions of reduced velocity; a) 
Imaginary part, b) Real part 


Increasing velocity produces a larger aerodynamic feedback. This aeroelastic coupling 
causes the dynamic mode frequency to decrease as velocity increases. This trend holds 
true until the system becomes unstable. This configuration destabilizes as a zero 
frequency root, aerodynamic in origin, migrates across the imaginary axis, that is, 
divergence occurs. 


The eigenvalues of the system for the divergence reduced velocity are distinguished in 
Figure 9 by solid squares. It is apparent that the dynamic mode still exists with a nonzero 
frequency when the system becomes unstable. At this velocity, the dynamic mode is a 
coupled structural and aerodynamic mode; the modal content and resultant system 
behavior will be addressed subsequently. 


Attention is now turned back to the aerodynamic roots, focusing on the real axis. The 
aerodynamic roots which lie on the real axis are of primary concern in the study of 
divergence. Two real poles originate from the present aerodynamic model. The 
existence of aerodynamic roots at zero velocity is addressed in the discussion portion of 
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this paper. For increasing airspeed, both roots initially become more stable. However, at 
approximately 75% of the divergence reduced velocity, one root changes direction and 
eventually destabilizes. The real part of the eigenvalues, shown in Figure 10, is 
indicative of the damping characteristics. The real part of a root that is losing damping 
will shrink. This indicates that energy is not being dissipated by this mode as effectively 
as at lower velocity. The real part of the dynamic root is becoming smaller, showing that 
it is losing its ability to dissipate energy. Once the static root destabilizes, the dynamic 
root, structural in origin, no longer tends towards instability. 


Modal Characteristics at Divergence: Non-dimensional Parametric Variations 

The configuration analyzed above exhibited divergence while the dynamic mode 
persisted. This behavior is predominant throughout much of the parametric space. This 
design space will now be explored. Variations in the structural non-dimensional 
parameters were performed. For each set of parameters, the aeroelastic equations of 
motion were constructed and eigenvalues found. At the divergence reduced velocity, the 
modal characteristics can be observed, specifically the frequency and damping of the 
mode which originated as the structural dynamic mode. 

A database of modal characteristics at divergence was generated to identify regions in the 
parameter space where the divergence mechanism changes from being associated with 
the structural root versus an aerodynamic root. The parameter variation results were also 
utilized in the design process for the wind tunnel configurations. The structural dynamic 
parameters varied in the database are elastic axis location, e/b, radius of gyration, r«, and 
mass ratio, p. 

The natural frequency of the pure structural torsional mode, 0l> (X , was also varied in the 
initial studies. Changing co u , however, was found to have no effect on the eigenvalue 
migration pattern as reduced velocity varied. From steady aerodynamic equations, it can 
be readily observed that to,* has no direct impact on the reduced velocity of divergence. 
Thus, the natural frequency is not studied in the parameter variation database. 

Varying the three parameters produces a four-dimensional parameter space. Thus, the 
results can not be shown by a single plot. Sample results are presented here as three- 
dimensional surface plots. For each surface, one of the parameters is held at a fixed 
value. Three surfaces are presented showing the ratio of the dynamic mode frequency at 
divergence normalized by the pure torsion frequency, C0o Q) (x . Three corresponding 
surfaces for the damping ratio at divergence, are also presented. 

The elastic axis position is Fixed for the surfaces shown in Figure 1 1 and Figure 12. 
Figure 1 1 presents a surface of frequency of the dynamic mode at divergence, normalized 
by the air off pitch frequency. The surface shows the variation of the divergence 
frequency as a simultaneous function of both radius of gyration and mass ratio. The 
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presented surface is for a nondimensional elastic axis position of 0.375, which 
corresponds to the wind tunnel model configurations. The solid square on this surface 
corresponds to wind tunnel configurations 2 and 3. The circle shown on the surface 
corresponds to wind tunnel model configuration #1, which has a higher radius of 
gyration. Note that the two wind tunnel model configurations lie in different plateaus of 
the surface. This difference will be shown to be indicative of a qualitative change in 
system characteristics. This figure is one slice from the parametric variation design 
space, here a four-dimensional space. As the elastic axis moves closer to the center of 
pressure, the surface becomes less smooth. The ridge that is shown in the back left 
corner of this parametric slice becomes a sudden hill. A trough develops in front of the 
hill and an additional ridge emerges which runs from low values of mass ratio diagonally 
across the space to low values of radius of gyration. The front right corner drops to a 
form a plateau where the frequency ratio becomes zero, or the dynamic mode has become 
real. 

Figure 14 presents the companion surface showing the damping ratio at the divergence 
condition. The damping information is presented as the angle between the imaginary axis 
and the eigenvalue. This angle is shown in radians, with a maximum magnitude of the 
angle of tc/ 2 or 1 .57 radians. 


0.6 



r 

u 


Figure 1 1 Surface of (Od/o^ as a function of mass ratio and radius of gyration. Elastic axis position is 
fixed, e/b=0.375. 
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Figure 12 Surface of as a function of mass ratio and radius of gyration. Elastic axis position is 
fixed, e/b=0.375. 


The radius of gyration is fixed for the surfaces shown in Figure 13 and Figure 14. Figure 
13 presents the surface of frequency of the dynamic mode at divergence, normalized by 
the air off pitch frequency for a radius of gyration of 0.5, which roughly corresponds to 
the primary wind tunnel model configuration. The surface shows the variation of the 
divergence frequency as a simultaneous function of mass ratio and elastic axis position. 
The solid square on this surface corresponds to wind tunnel configurations 2 and 3. 

Again, this figure is one slice from the parametric variation design cube. As the radius of 
gyration decreases, the surface becomes overall smoother. As the radius of gyration 
increases, the trough in the middle of the mass ratio range migrates to lower values and a 
new trough forms at high values. The new trough at higher mass ratios grows as radius 
of gyration increases and the overall surface resembles two plateaus connected by a steep 
grade. The plateau in the frequency ratio surface at high mass ratios has a value of zero, 
where the dynamic mode has become real. Figure 14 presents the companion surface 
showing the damping ratio at the divergence condition. 
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Figure 13 Surface of oW^ as a function of mass ratio and elastic axis position. Radius of gyration is 
fixed, r„=0.5. 
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Figure 14 Surface of ^ as a function of mass ratio and elastic axis position. Radius of gyration is 
fixed, r u =0.5 

The mass ratio is fixed for the surfaces shown in Figure 15 and Figure 16. Figure 15 
presents the surface of frequency of the dynamic mode at divergence, normalized by the 
air off pitch frequency for a mass ratio of 50, which roughly corresponds to the wind 
tunnel model configuration 2. The surface shows the variation of the divergence 
frequency as a simultaneous function of elastic axis position and radius of gyration. The 
solid square on this surface corresponds to the wind tunnel configurations 2 and 3. The 
rolling hills at low values of elastic axis location and the gentle slope into a single plateau 
at high values of elastic axis location are characteristic throughout the range of mass ratio 
considered, which was from 20 to 200. Lower values of mass ratio make the hills more 
dramatic, while higher mass ratios make the valleys sink to zero. Figure 16 presents the 
companion surface showing the damping ratio at the divergence condition. 
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Figure 16 Surface of as a function of elastic axis position and radius of gyration. Mass ratio is 
fixed, p=50. 


As seen in the sample parameter surfaces, there are regions where the structural dynamic 
roots have become real, indicated by the frequency ratio becoming zero. These regions 
indicate where the divergence mechanism will look like the traditional interpretation of 
divergent behavior. This information occurs in a relatively small region of the 
nondimensional parameter space, however. This region is shown in Figure 17. 
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Figure 17 Non-dimensional parameter space where traditional divergence mechanism 
occurs 
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From the parametric database, many observations can be made. The surfaces of 
frequency ratio and damping at divergence are neither uniform nor monotonic. However, 
in general locating the elastic axis just aft of the center of pressure tends to make the 
typical section diverge in a traditional manner. Large radii of gyration also produce this 
effect, as do large mass ratios. 

The three-dimensional parameter variations reveal that the parameter space is divided 
into three distinct regions. 1 ) The structural dynamic root migrates to the real axis as the 
reduced velocity increases. This mode diverges, as in the traditional interpretation of 
divergence. 2) A real aerodynamic eigenvalue diverges. The structural dynamic mode 
still exists as a complex mode at the corresponding reduced velocity. 3) A real 
aerodynamic root diverges. The structural dynamic root has previously become real, 
migrated further left along the real axis, becoming more damped, and then becoming 
complex again prior to divergence. 

The wind tunnel model configurations were designed to fall into the second category and 
demonstrate divergence of a static mode whose origin lies in an aerodynamic root. 

Three wind tunnel model configurations will be discussed in this paper. Analytical 
stability results for what will come to be known as configuration # 2 have already been 
presented. Each of these configurations diverges as a static mode which originated in the 
aerodynamic model becomes unstable, while the dynamic mode with its origin in the 
structural modelpersists at a non-zero frequency. The structural dynamic parameters for 
all three configurations are provided in Table 5. For definitions of non-dimensional 
parameters, see Table 3. 


Config 


K„ 

C0„ 

fa 


r« 


# 

(slinch- 

. 'I 

in') 

(lbr 

in/deg) 

(rads/sec) 

(Hz) 




1 

0.1147 

.90 

21.2 

3.37 

0.0046 

0.741 

107.9 

2 

0.021 

0.90 

49.5 

7.88 

0.0053 

0.459 

51.4 

3 

0.021 

2.78 

87.1 

13.86 

0.0035 

0.462 

50.8 


Table 5 Structural dynamic parameters associated with wind tunnel model configurations 


Configuration #1 has the same torsional stiffness as configuration #2, but the trailing 
edge segment is made of Tungsten, which substantially increases the pitch inertia. The 
eigenvalue migration is qualitatively different than that presented for configuration #2. 
The root locus will be presented shortly. Configuration #3 has the same pitch inertia as 
configuration #2, but the torsional stiffness was increased. The eigenvalue pattern is 
identical to that of configuration #2. The change in torsional stiffness manifests itself in 
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changing the divergence dynamic pressure of the system, but the non-dimensional locus 
does not change at all. 

Table 6 lists the analytical calculations for divergence conditions for the three 
configurations. Note that the non-dimensional divergence condition is invariant with 
torsional stiffness. This is shown by the agreement of the reduced velocity for 
configurations 2 and 3. Also note that the physical divergence condition is invariant with 
pitch inertia. This is shown by the agreement of the velocity or dynamic pressure for 
configurations 1 and 2. 


Configuration # 

Reduced 

Velocity 

Velocity 

Dynamic Pressure 



(in/sec) 

(mph) 

(psf) 

(N/nr) 

1 

8.89 

754 

42.8 

4.6 

222 

2 

3.8 

754 

42.8 

4.6 

222 

3 

3.8 

1324 

75 

14.25 

687 


Table 6 Analytical calculation of divergence conditions 


The divergence conditions calculated using the aeroelastic eigenanalysis can be compared 
to the divergence conditions calculated using the equations of static equilibrium of the 
system. The equations of static equilibrium for the single degree of freedom typical 
section can be written by equating the aerodynamic and structural moments which act at 
the elastic axis. Equation 1 1 . The aerodynamic moment for a symmetric airfoil can be 
expressed in terms of the lift curve slope and total angle of attack, which is comprised of 
the rigid angle of attack, Oku and the elastic increment, Ok, Equation 12. The structural 
restorative moment. Equation 13, is proportional to the elastic increment. Setting them 
equal and rearranging produces a ratio of elastic increment to rigid angle of attack. 
Equation 14. For a finite rigid angle of attack, the elastic increment will become infinite, 
diverge, if the denominator of the right hand side becomes zero. This provides the 
expression in Equation 15 for calculating the divergence dynamic pressure. Recasting 
the equation in non-dimensional quantities produces Equation 16. 

M 4 = Mg Equation 11 


M a = (0() + (Z e ) Equation 12 


M s =K a a e Equation 13 
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Equation 14 


a e _ <I SeC L a 
«() K a -qSeC L) 


Qd ~ 


K 


a 


SeC I 


Equation 15 





Equation 16 


The divergence conditions were calculated lor the three wind tunnel model 
configurations, utilizing the parameters listed in Table 5. The results, summarized in 
Table 7, compare almost exactly with the analytical results shown in Table 6. 


Configuration 

Divergence Dynamic 
Pressure 

(psf) 

Divergence Reduced 
Velocity 

1 

4.69 

8.89 

2 

4.69 

3.80 

3 

14.49 

3.80 


Table 7 Static equilibrium calculations of divergence conditions 


Divergence of an eigenvalue which originates in the aerodynamic model is contrary to 
the traditional interpretation of system behavior at divergence, although several similar 
phenomena have been reported by earlier researchers, as discussed in the introduction 
and background sections of this thesis. The modal content of the instability and the 
dynamic mode are addressed in the following study of the eigenvectors which provides a 
deeper insight into the observed phenomena.. 


Eigenvector Study 

The eigenvectors of the aerodynamic and/or aeroelastic system are now examined in 
detail. While it is common to examine the eigenvalues for information on system 
behavior and stability characteristics, it is quite uncommon to attempt to garner insights 
from the eigenvectors. As noted by Bisplinghoff and Ashley': “(the aeroelastic 
eigenvalues) have associated with them eigenfunctions when the complex representation 
is used. Since (the absolute magnitude of eigenfunctions, or eigenvectors) may be 
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specified arbitrarily, the mode shape is completely defined by the amplitude ratio and a 
phase angle between . . . degrees of freedom. This information has only minor interest in 
stability studies, and modes are often not even calculated. Their variation with airspeed 
can throw some light on the physical nature of (the instability), however.” 

Determination of the flutter mode shape is the most common use of the eigenvectors 
resulting from an elastic analysis. The eigenvectors contain a wealth of information that 
is only hinted at by the eigenvalues. An individual eigenvector can reveal the frequency 
and damping of the associated eigenvalue. It also provides a ratio of the energy present 
in the mode due to each component of the state vector. In the case of the eigenmodes 
analyzed here, an eigenvector provides information on the modal content: whether a 
mode should be considered as primarily structural, primarily aerodynamic or as a hybrid, 
aeroelastic mode. The orthogonality or lack of orthogonality among the eigenmodes 
indicates whether energy can be transferred from one mode to another. This can be 
important in understanding when, how and why an aeroelastic system destabilizes. 


Aerodynamic Eigenmodes 

The eigenmodes of the aerodynamic system are considered first. These modes do not 
change as velocity increases, but it is instructive to know what the modes look like. The 
aerodynamic eigenmodes contain the modal vorticity for each aerodynamic element. 
Eight of the eigenmodes are shown in Figure 1 8. The eigenmodes are presented for an 
aerodynamic model with 10 elements on the wing and 180 elements in the wake. The 
modal vorticities are plotted at the chord-wise location of the associated aerodynamic 
elements. 

The first mode, a real mode, resembles a static pressure coefficient distribution over the 
wing, with little participation from the wake. The second mode is also a real mode 
resembling a static pressure distribution over the wing. However, this mode contains a 
large amount of wake participation. The remaining aerodynamic modes are complex and 
are comprised primarily of an oscillating wake. The wing vorticities are insignificant 
compared to those in the wake. The modes are ordered by increasing frequency. Each 
mode contains a single frequency; as the frequency increases or as the mode number is 
advanced more oscillations are observed in the wake. 
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Figure 18 Selected aerodynamic eigenmodes 


Aeroelastic Eigenvectors 

The aeroelastic eigenvectors are now studied from several perspectives. The first 
approach taken in examining the aeroelastic eigenvectors is to study the behavior 
associated with individual modes. Following this, the relationship between two 
eigenvectors is examined. 

In a numerically stiff set of ordinary differential equations 4 , the system behavior is seen 
to be dominated by the lightly damped and unstable modes. The disparity in the time 
scales of components of the system allows the overall behavior to be studied by 


4 Kincaid. D.R.. and E. W. Cheney, Numerical Analysis: Mathematics of Scientific Computing 
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observing only a few eigenmodes of the system. Thus, in this discussion of system 
behavior, the modal participation factors associated with only a few modes are examined. 
Three eigenvectors are considered: 1 ) the dynamic mode which originates as the 
structural dynamic mode; 2) the least stable real mode which originates in the 
aerodynamics and becomes the source of divergence and 3) the second real aerodynamic 
mode. This information is repeated for reference in Table 8. 


Eigenvector Numbering 

Real or Complex 

Origination 

1 

Complex 

Structural dynamic mode 

2 

Real 

Aerodynamic mode 

3 

Real 

Aerodynamic mode 


Table 8 Eigenvector numbering and description 


Examination of individual eigenvectors > 

The first approach taken in examining the aeroelastic eigenvectors is to study the 
behavior associated with individual modes. The eigenvector associated with a particular 
eigenvalue can be viewed as the set of modal participation factors for each degree of 
freedom. Note that the eigenvectors are invariant under the transformation from discrete 
to continuous time domain. A proof of this is given in Appendix B. The dynamic mode 
and the destabilizing static mode are examined in detail; the vorticity portion of the 
eigenvectors is emphasized. 

The analysis results presented here are for wind tunnel model configuration #2. The 
eigenvectors have been normalized to have unity magnitude and phased such that the 
structural dynamic generalized displacement coordinate, a, has zero phase. 

The dynamic mode near zero velocity is considered first. The modal participation at a 
low reduced velocity, V=0.225, is presented for the dynamic mode in Figure 19. At this 
velocity, the mode is almost a pure structural pitch mode. The associated eigenvalue is 
identified in the continuous time root locus, Figure 9, by the diamond symbol. The real 
and imaginary parts of the modal participation are plotted as a functions of chord-wise or 
downstream position. At this low velocity, the aerodynamics are being driven at the 
frequency of the structural mode. The portion of the eigenvector associated with the 
vorticity at each aerodynamic control point, referred to as the vorticity participation, 
shows that most of the aerodynamic energy associated with this mode is in the wake. The 
first ten participation factors correspond to elements on the airfoil. Only these vorticities 
can produce forces on the airfoil. At this velocity, there is very little aerodynamic energy 
being imparted to the airfoil. 
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Normalized Chordwise Location, (x/b ) 

Figure 19 Eigenvector associated with dynamic mode near zero velocity, associated continuous time 
eigenvalue X = -.16 + j 49.2 

The wake portion of the vorticity participation appears as a negatively damped sinusoid 
when viewed spatially. The eigenvector provides a snapshot of the vorticity distribution. 
Initial examination of the data may lead one to conclude that the system is unstable. In 
fact, the opposite is indicated. For a stable system, the vorticity being shed from the wing 
into the wake will decrease as time advances. The vorticity on the last wake element at 
time n is the same as the vorticity on the first wake element at time n-N wa k C . Thus, the 
spatial vorticity distribution could also be thought of as a time history, where time 
originates at the wake trailing edge and proceeds towards the airfoil. 

Near the divergence reduced velocity, the eigenvector associated with the dynamic mode 
contains significant participation from both the structural dynamic and the aerodynamic 
states. Figure 20 shows the vorticity participation spatially for a velocity just above 
divergence, V= 3.85. The number of oscillations to be expected in the wake, N cyc i es , can 
be estimated using the frequency of the associated eigenvalue, 0) mo de, the reduced 
velocity, V, and the number of aerodynamic discrete elements in the wake and on the 
airfoil, N wa k C and M: 

\r _ ^mode ^ wake 

” cycles ~ TTT , — Equation 1 7 

J (0 0 MkV 

Using the values for the divergence condition results in a prediction of 0.48 spatial 
oscillations; the vorticity participation in Figure 20 is consistent with this estimate. 
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Figure 20 Eigenvector associated with dynamic mode at the divergence velocity, associated 
continuous time eigenvalue A. = -17.7 + j 25.9 

It is also instructive to view the eigenvectors in terms of the magnitude and phase angle. 
The vorticity portion of the dynamic mode eigenvector is presented in this format in 
Figure 21. For 5 values of reduced velocity, listed to the left of each magnitude plot, the 
magnitude and phase of the eigenvector components are plotted as a function of the 
chord-wise location of the aerodynamic box. The modal vorticities on the wing are 
shown by the circles. The wake modal vorticities are shown by dots, which appear as 
solid lines due to the dense spacing. The last element of the wake has not been shown in 
these plots- it will be discussed separately. 

The magnitude plots show that as the velocity is increased and approaches the divergence 
speed, more modal energy is contained in the aerodynamic portion of the eigenvector. 
The divergence reduced velocity for this configuration was calculated as 3.8. The 
magnitude plots indicate that beyond this velocity, the wing vorticity participation 
decreases. As previously discussed, the vorticities on the wing determine the importance 
of the aerodynamic feedback. As velocity increases, the aeroelastic coupling in the 
dynamic mode increases as evidenced by the growing magnitudes of wing vorticities. 

The phase plots also provide much useful information. The dynamic mode eigenvector 
can be viewed as if the aerodynamics are being forced at the modal frequency. As the 
airspeed advances, the frequency of the excitation changes. As shown in the cases of the 
zero airspeed and divergence, the number of cycles expected in the wake can be 
approximated using Equation 17. The phase information from Figure 21 is summarized 
in Table 9. The cycles in the wake are estimated and an approximate value of the 
frequency is calculated. These are shown in the table to compare almost exactly with the 
imaginary parts of the associated eigenvalues. 
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Figure 21 Vorticity portion of dynamic mode eigenvector for several velocities, magnitude and phase 
as functions of chord-wise location 


Reduced Velocity 

Number of cycles 
counted in the 
wake 

(Omode calculated 
from Equation 17 
(rads/sec) 

(Omodc from 
aeroelastic 
eigenvalue 

(rads/sec) 

1 

5.5 

48 

48.1 

2 

2.6 

45 

44.5 

3 

1.4 

36 

36.7 

3.5 

1.0 

30 

29.9 

5 

.56 

24 

24.1 


Table 9 Dynamic mode frequencies estimated from wake portion of dynamic mode eigenvector and 
calculated from analysis 
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The aeroelastic system studied destabilizes as a real eigenvalue moves into the right half 
plane. The vorticity participation factor associated with this mode resembles a pressure 
coefficient distribution on the airfoil elements, while the wake contains almost no 
participation except tor the last element. The vorticity participation factor at an example 
reduced velocity, chosen here to correspond with divergence, V= 3.85, is presented in 
Figure 22. As the reduced velocity changes, it is the participation of the last wake 
element is especially interesting. The magnitude and phase of this element of the 
eigenvector is plotted versus reduced velocity in Figure 23. Note that these eigenvectors 
have an overall magnitude of 1 . Initially, nearly all of the vorticity participation resides 
in the last element of the wake. As velocity increases, all of the wake elements begin to 
participate in the mode. Just prior to divergence, the participation of the last wake 
element drops sharply. At the divergence velocity, all of the vorticity participation is on 
the airfoil; the wake factors are zero. As the system moves beyond the divergence 
velocity, the behavior of all of the vorticity participation factors change. The last wake 
element quickly becomes influential again, but now with vorticity that is negative, or out 
of phase, with the airfoil vorticity. As velocity is further increased, the participation of 
the last wake element smoothly, asymptotically, approaches zero. Also beyond 
divergence the overall wake vorticity participates. 



Normalizad Chordwise Location, (x/b) 


Figure 22 Eigenvector associated with the unstable static mode just above the divergence velocity 

Transition from stability to instability produces dramatic changes in the associated 
eigenvector. While the eigenvalue smoothly traverses across the imaginary axis, the 
character of the vorticity participation changes sharply. 
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Figure 23 Vorticity participation of last wake element, associated with the real eigenvalue 
that destabilizes 



Vorticity Ratios as Indicators of Modal Participation 

Each of the modes of an aeroelastic system contains structural and aerodynamic 
participation to varying degrees. This participation changes as the velocity increases and 
the feedback from the aerodynamics into the structure increases. Furthermore, the 
structural contributions and the aerodynamic contributions are different for different 
configurations. As an indicator of the aerodynamic contribution to the dynamic mode 
and the divergence mode, the modal vorticities on the wing were summed. These sums 
were normalized by the summation of the wing and wake modal vorticities. This ratio 
indicates the amount of aerodynamic participation in the mode. 

Four configurations are examined in Figure 24. The root loci are shown in the left 
column. The roots are shown which correspond to the structural dynamic mode, which 
originates as a complex pair, and also the two real aerodynamic eigenvalues. The zero 
airspeed values for each of the roots are shown in the plot with open triangles. As the 
reduced velocity is increased towards the divergence condition, the roots of the system 
are shown by the dots. The value of each root at the instability dynamic pressure is 
denoted by a solid triangle. The vorticity ratios are shown in the right column. The 
vorticity ratios are presented for the dynamic mode which is structural in origin, (mode 
I ), and the divergent static mode which is aerodynamic in origin, (mode 2), as functions 
of reduced velocity. The area above each curve represents the energy dissipated into the 
wake. The area below each curve represents the modal energy imparted to the structure. 
A higher curve indicates that a larger portion of the modal energy can be attributed to the 
aerodynamics. 
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Figure 24 Root loci and vorticity ratio plots for a) wind tunnel model configuration 2; b) wind tunnel 
model configuration 1; c) the intermediate configuration and ; d) the traditional divergence 
configuration 
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The first configuration is the previously analyzed wind tunnel model configuration #2. 
The dynamic mode persists for this configuration at divergence. The root locus shows 
the dynamic mode frequency at greater than half the natural frequency of the air off 
system. This suggests that the structural influence is a dominant participant in the mode. 
The vorticity ratio tells a similar story. The dynamic mode’s vorticity ratio is shown by 
the smaller symbols. The aerodynamic contribution to the mode is low throughout the 
reduced velocity range. Thus, the dynamic mode is primarily a structural dynamic mode. 
Revisiting the root locus, it can be observed that the primary contribution of the 
aerodynamics to this mode is a large amount of damping. The vorticity ratio associated 
with the static divergence mode is also shown. At divergence, the aerodynamic 
participation is shown to increase dramatically and dominate the mode. Consideration of 
the root locus and the vorticity ratio indicates that divergence for this configuration is 
primarily aerodynamic. 

The second configuration shown in Figure 24 corresponds to wind tunnel model 
configuration #1 . The frequency of the dynamic eigenvalue has decreased significantly 
at divergence compared to the previous configuration. The vorticity ratio indicates that 
the dynamic mode has more aerodynamic participation than the previous configuration. 
The dynamic mode is still observed to have a nonzero frequency at divergence and be 
primarily driven by the structural participation. The divergent mode is primarily 
aerodynamic in nature. 

The third configuration shown does not correlate with a constructed wind tunnel model 
configuration. It is analyzed to show' the progression of the aerodynamic participation in 
the dynamic mode. The vorticity ratio indicates that the dynamic mode has more 
aerodynamic participation than the previous configurations. 

The fourth configuration illustrates the divergence mechanism that is traditionally 
envisioned. The parameters used in this example do not represent a buildable 
configuration with known materials, but are presented to show that the analysis 
methodology is not single-minded. This configuration was generated by modifying the 
parameters associated with wind tunnel model configuration #1. The elastic axis was 
moved to X A inch aft of the center of pressure. The radius of gyration was then doubled. 
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Description 

Parameter 

Value 

Units 

Semi-chord 

b 

4.0000 

inches 

Span 

span 

21.0000 

inches 

Radius of gyration 


1.6034 


mass ratiio 

n 

108.0000 


torsion mode frequency 

co« 

19.6000 

radians/second 

Elastic axis location 

e/b 

0.0561 


Number of aerodynamic 
elements on the wing 

M 

10.0000 


Total aerodynamic 
elements 

N 

180.0000 


density of air 

Pair 

0.0000 

slinches/inch' 

aerodynamic relaxation 
factor 

a 

0.9960 


size of aerodynamc 
element 

Ax 

0.8000 

inches 


Table 10 Parameter Values for traditionally divergent configuration 

The corresponding root locus, which is presented in the last row of Figure 24, indicates 
that the structural dynamic originated mode frequency decreases with increasing reduced 
velocity until they turn from a complex pair of roots into two real roots. Shortly after the 
pair become real, one root becomes more highly damped and one becomes unstable. As 
the reduced velocity is increased, the real aerodynamic eigenvalue migrates toward the 
left, or damped condition. The structural dynamic eigenvalues migrate towards the real 
axis. As they hit the axis, the complex pair of eigenvalues both become real. One of 
them becomes more highly damped and the other one becomes unstable. 

The non-dimensional parameters used in this last example case can be produced by 
physical parameters which approach buildability only for semichord values which exceed 
wind tunnel blockage guidelines. Also, there are a few additional issues regarding the 
potential for successful fabrication and testing of this configuration. In locating the 
elastic axis so near the center of pressure, there is very little room for error or change in 
the fabrication process. The dynamic pressure of the instability is proportional to the 
inverse of the distance from the center of pressure to the elastic axis. If a 1/16 inch error 
exists in the location of the pivot axis, the divergence dynamic pressure would be 75% of 
the anticipated value. A second consideration is the potential for a camber mode to be 
induced. The airfoil is a 1 /32-inch thick aluminum shell. At the trailing edge, a tungsten 
mass is attached with set screws. Moving the rotational axis forward, which is required 
to generate the traditional mechanism configuration, means moving the rotational axis 
away from the supporting, stiffening interior spars and also giving the tungsten a longer 
moment arm. These factors will make the airfoil section tend to deform in the 
streamwise direction. All analyses to date have been performed assuming a rigid airfoil. 
The airfoil section, however, is a closed cell. This means that it is a fairly stiff structure, 
and the deformation may not be a cause for concern. 
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Modal Moment Comparison 


The structural and aerodynamic moments can also be calculated by using the eigenvector 
information. The static structural moments corresponding to each mode were calculated 
using the static moment equation. Equation 1 3, and the modal results for angular 
displacement. The static aerodynamic modal moments are calculated in similar fashion, 
employing the static aerodynamic equation. Equation 12. The results are compared in 
Figure 25. In this comparison, the angle of attack component of the dynamically 
determined eigenvectors was employed. For each velocity, the ratio of the aerodynamic 
to static moment is a ratio of dynamic pressure to divergence dynamic pressure. This can 
be shown by comparison of Equation 2 and Equation 3. Subcritically, for each mode, 
the static structural modal moment is less than the aerodynamic static moment. Thus, 
each mode has enough structural restorative power to counteract the effect of the 
aerodynamic moment, and the modes are stable. In the supercritical case, the static 
mode, labeled mode 2, indicates that the aerodynamic moment is too large for the 
structure to restore the system to equilibrium. This indicates that the mode is statically 
unstable, as indicated previously by the eigenvalue analysis. 
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Figure 25 Modal moments, configuration #2 


V 

V 

V 

V 

V 


V 

V 


V 

V 


1 1 V* 


fj 



53 




The dynamic modal moments associated with any eigenvector must have equal structural 
and aerodynamic components as defined in the dynamic equations. This modal quantity 
is also shown in the figure for both of the modes represented. A difference between the 
static and dynamic moment is the inclusion of the oscillatory portion of the motion. For 
the dynamic mode, particularly subcritically, there is a large oscillatory component. This 
is illustrated by the large different between the static structural moment for mode 1 and 
the dynamic moment for the same mode. Another significant difference between the 
static moment results and the dynamic moment results occurs supercritically. The static 
equations essentially enforce a neutral stability assumption. Because the static mode, 
labeled mode 2. is unstable, there is a large difference between the static and dynamic 
moments. 


Orthogonality between eigenvectors 


Comparisons were made between pairs of eigenvectors, employing the techniques 
demonstrated by Afolabi, Pidaparti and Yang 5 . 

The eigenvectors identified in Table 8 were compared, two at a time by finding the phase 
angle between them. Modes which are in phase, or have 0° separating their orientations, 
will tend to feed energy into each other and potentially amplify the motion. Modes which 
are out of phase, or have angles of 180°, will tend to act against each other, canceling out 
the energy and motion of each other. Orthogonal modes, generally thought to be 
incapable of exchanging energy from one mode to another, would be at 90° to each other. 
Near a modal coalescence, a loss of orthogonality could be expected to occur. The angle, 
0, between the eigenvectors, ij/i and \j/„ is computed as the inverse cosine of the 
normalized inner product. 


&ij =cos 


Vi ' 

JMM, 


Equation IS 


For complex modes, inner product is seldom real, resulting in an inability to compute the 
arccosine. This difficulty is avoided by employing “realification” of the complex 
eigenvectors. Realification of a complex array is a stacking of the real parts and then the 
complex parts, turning a vector of length n into a vector of length 2n. For either the real 
or the realified complex case, this procedure is essentially finding a weighted average of 
the phase of all eigenvector components. 


In analyzing the current aeroelastic system, the coordinates are not all of comparable 
physical or mathematical quantities- some being structural dynamic generalized 
coordinates and some being aerodynamic vorticities. This makes interpretation of the 


5 Afolabi, Dare, Ramana M. V. Pidaparti, and Henry Y. T. Yang, Flutter Prediction Using an Eigenvector 
Orientation Approach 


54 



results more subtle and difficult. Figure 26 shows the phase angles between the 
eigenvectors as functions of reduced velocity. All components of the eigenvector were 
utilized in this comparison. 



Figure 26 Angle between eigenvectors as a function of reduced velocity, using all elements of the 
eigenvector 

The angle between the two real modes’ eigenvectors, denoted 2 and 3, is shown in the 
figure with a dashed line connecting circular symbols. At low velocities, the 
eigenvectors are nearly out of phase. The angle between the two real modes shows a 
sudden transition initiated just prior to divergence. At the divergence velocity, the two 
real eigenvectors are orthogonal; above the divergence velocity, the phase overshoots 90° 
and then asymptotically reapproaches 90°. 

The angles between the real modes and the dynamic mode are oscillatory as reduced 
velocity increases. The angle between the dynamic mode and the mode, which 
destabilizes is shown in the figure with a solid line connecting square symbols. The 
angle between the dynamic mode and the stable aerodynamic-originated mode is 
indicated by a dotted line connecting triangular symbols. For both comparisons, the 
frequency of the oscillations decreases and the magnitude increases. The underlying 
modeling which produces these characteristics is addressed by separately considering the 
eigenvector in portions corresponding to the wake vorticity components, the wing 
vorticity components and the structural dynamic generalized coordinates. 

The vorticity in the wake is considered first. Details of the vorticity portions of the 
dynamic mode eigenvector were previously presented for several reduced velocities. 
Figure 19, Figure 20 and Figure 21. The angles between the eigenvectors were computed 
as previously, except that only the subset of the eigenvector components corresponding to 
vorticities in the wake were included. As before, the last wake element has been ignored. 
Figure 27 shows the angles between subsets of the eigenvectors as functions of reduced 
velocity. There are oscillatory patterns shown relating the phase angle of the dynamic 
mode eigenvector to both real eigenvectors. As in the case which utilized all components 
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of the eigenvectors, the frequency decreases and the magnitude increases as reduced 
velocity increases. From the detailed eigenvector plots of the dynamic mode, it can be 
seen that the cycling of the phase is more rapid at the lower velocities. The growing 
magnitude indicates the same thing. The more cycles that the phase goes through, the 
smaller the average phase of the vector. For the comparison with the stable real 
eigenvector (eigenvector 3), the phase angle smoothly oscillates for the entire range of 
velocity presented. The phase angle between the dynamic mode and the unstable mode, 
however, undergoes a change near the divergence velocity. The slope of the phase angle 
curve changes sign. This is attributed to the change in sign of the real mode eigenvector. 



Figure 27 Angle between eigenvectors as a function of reduced velocity, using wake vorticities 
without last element 

Lack of orthogonality serves to indicate the potential for energy being passed from one 
mode into another. The wake portion of the eigenvector does not have a direct impact on 
the aeroelastic feedback forces generated. To offer insight into this behavior, the wing 
vorticity and structural dynamic portions of the eigenvectors are examined. Because of 
the mismatch in physical quantities represented, they are considered separately. 

The vorticity on the wing is examined next. The angles between the eigenvectors were 
computed, as before, except that only the portions of the eigenvectors corresponding to 
vorticities on the wing were included. Figure 28 shows the angles between subsets of the 
eigenvectors as functions of reduced velocity. The real modes are seen to be in phase 
throughout the velocity range presented. The relationships of the dynamic mode to the 
real modes are nearly identical, as would be expected after examining the phasing 
between the real modes. Near zero velocity, the dynamic mode is nearly orthogonal to 
the real modes. This orthogonality is quickly lost as the airspeed is increased. This 
indicates that the modes can not exchange energy when there is no velocity. From a 
physical standpoint, this indicates that there is no aeroelastic feedback, or interaction of 
the structural and aerodynamic entities, when the velocity is low. As the airspeed 
increases, the angle stays relatively constant. The upslope that starts just prior to 
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divergence is attributed to redistribution of vorticity between the real modes, which is 
also indicated by examining the angle between eigenvectors 2 and 3. 



Figure 28 Angle between eigenvectors as a function of reduced velocity, using wing vorticities only 


The structural dynamic generalized coordinate contributions to the eigenvectors are 
examined next. The angles between the eigenvectors were computed, using only the 
portions of the eigenvectors corresponding the structural portion. Figure 29 shows the 
angles between subsets of the eigenvectors as functions of reduced velocity. It has been 
mentioned that the eigenvectors were normalized such that the angle of attack generalized 
coordinate has a phase of 0°. This is true for each eigenvector individually examined. 

The structural dynamic eigenvector segment consists of this angular displacement and the 
velocity of this coordinate. For this reason, the angle between the eigenvectors, which is 
a weighted average for all components of the eigenvector included in the analysis, does 
not have a phase of Of. The angle between the real modes, indicated by the dashed line 
connecting circular symbols, begins at low velocity with the eigenvectors nearly in phase. 
At divergence, there is a sharp transition, which shows that they are orthogonal at 
divergence, and are out of phase for velocities above divergence. The orthogonality at 
divergence indicates that modal energy can not be transferred from the unstable mode 
and dissipated by the stable real mode. For velocities below the divergence speed, the 
modes are nearly orthogonal, indicating that little energy is transferred between modes 
through the structural dynamic participation. As the velocity approaches divergence, the 
angle between the dynamic mode and the unstable mode changes. The modes start to 
lose their orthogonality. Just prior to divergence, the separation angle is approximately 
70°, so some structural dynamic energy can be exchanged between the modes. At 
divergence, the modes are orthogonal, so no energy is exchanged between them. Beyond 
divergence, the modes are tending towards being out of phase with each other, indicating 
that the modal motions would oppose each other. Meanwhile, the stable real 
aerodynamic mode and the dynamic mode are becoming more in phase. Coupling 
between these two modes is highly likely at higher velocities. This coupling can be 


57 




observed also from the root locus of the eigenvalues for velocities above divergence. The 
migration pattern of the dynamic mode eigenvalue seen in Figure 9 is clearly influenced 
by the presence of the stable real mode. 



Figure 29 Angle between eigenvectors as a function of reduced velocity, only the structural dynamic 
portion of eigenvectors used 

The eigenvectors were scaled such that the angular displacement has a phase ol zero. It 
makes sense then, to examine the structural dynamic portion of the eigenvectors by 
considering only the contribution of the angular velocity. These results are shown in 
Figure 30. The phase indicates that, at divergence, the unstable mode’s eigenvector 
phase changes by 180°. The dynamic mode starts orthogonal to the real modes, at 90°. 
This indicates that the dynamic mode either leads or lags the real modes of the system, 
such that the modes are orthogonal. As the reduced velocity increases, the phase between 
the dynamic mode and both real modes tend towards becoming in phase. At divergence, 
however, the curves separate. The angle between the structural dynamic mode and the 
unstable mode changes suddenly by 90°. The slope of the curve also changes such that it 
is now increasing with increasing velocity. As velocity increases, the increasing phase 
difference indicates that the modes tend towards becoming out of phase; they will no 
longer interact to accentuate each other’s motion. The phase difference between the 
stable mode and the dynamic mode continues its migration towards 0°, indicating that the 
modes are becoming more in-phase and the will be more apt to exchange modal energy at 
post-divergence velocities. 
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Reduced Velocity 

Figure 30 Angle between eigenvectors as a function of reduced velocity, only angular velocity 
component 



Figure 31 Angle between eigenvectors as a function of reduced velocity, using only last element of the 
wake 
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Computational Issues for Simultaneous Solution of Aerodynamic and Structural 

Equations 


Transformation Compatibility 


To incorporate the discrete time aerodynamic model into aeroelastic equations, the 
structural dynamic model must be cast in discrete time also. The structural dynamic 
equations contain first and second derivatives that could be approximated using a central 
difference technique. While this is convenient and easy, this method results in a 
mismatch of discrete time transformations. Central differencing produces discrete time 
equations to which a first order Tustin transformation' , Equation 19, 
must be applied to obtain the proper continuous time results. 

X — — — — — ■ Equation 19 

At (Z + l) 


The Tustin transformation is equivalent to the first term in a series expansion of the zero 
order hold transformation presented in equation 2. In these transformations, the sample 
interval, At, establishes the relationship between the discrete time eigenvalues, z, and the 
continuous time eigenvalues. A,. The aerodynamic equations which were generated with a 
zero order hold discretization, are solved simultaneously with the discretized structural 
dynamic equations. Thus, it is desirable to have structural dynamic equations that would 
also be correct when a zero order hold transformation is applied. This is easily 
accomplished through standard discretization techniques 6 . Accepting the mismatch in the 
transformations results in a phenomenon that resembles aliasing. However, as the time 
step becomes small, the zero order hold transform and the Tustin transform become 
approximately equivalent. 


Aliasing 

The equations have been constructed in the discrete time domain. Given data at discrete 
times, a transformation can be utilized to approximate the response in continuous time. 
There are limitations to discrete time transformation methods, aliasing is the primary 
concern 3 7 To avoid aliasing, a continuous time signal must have 2 samples per period of 
period of the highest frequency to be resolved. The aerodynamic equations arose from 
the fundamental concept of vorticity being convected downstream at a velocity, U. The 
equations are valid only if the relationship U=Ax/At is maintained. It is thus observed 


6 Phillips, Charles L„ and H. Troy Nagle, Jr, Digital Control System Analysis and Design. 
1 Oppenheim, Alan V., and Ronald W. Schafer, Discrete-time Signal Processing 
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that the minimum velocity, at which the system may be accurately analyzed, is set by the 
spatial discretization and the maximum frequency that is important to the problem. 
Another interpretation is that for frequency and velocity ranges of interest, the minimum 
number of aerodynamic elements required to avoid aliasing can be approximated. This 
can serve as a guideline in selecting the spatial discretization required for a given 
problem. There are additional implications of the discrete time effects when the 
aerodynamic equations are combined with the structural dynamic equations or control 
laws. 


Methods of Stability Analysis 

The aeroelastic stability analyses, which require variation of the velocity, were performed 
using a single spatial aerodynamic discretization. This was accomplished by adjusting 
the temporal discretization to produce the proper velocities. There are several 
complications in performing the analyses in this manner: ( 1 ) a separate transformation 
rule must be applied for each velocity; and (2) interpreting the discrete time eigenvalues 
is not intuitive. The aerodynamic matrices are unchanging for different velocities, but the 
matrices which couple them to the structural dynamics are not. The resulting aeroelastic 
eigenvalues change with each velocity. The migration of the eigenvalues in the discrete 
time domain is not due solely to the velocity change, but to a combination of velocity and 
sample rate change. 

A brief study was conducted to look at the results when a consistent sample rate was 
utilized, meaning that as the velocity changed, the spatial discretization changed. This 
required constructing a new aerodynamic model at each velocity. There was negligible 
effect on the continuous time eigenvalues. The discrete time eigenvalue pattern 
associated with the structural dynamic mode changed significantly. It was observed, 
however, that the discrete time eigenvalue pattern in this case is nearly identical to the 
pattern produced when the eigenvalues from the nominal analysis method are 
rediscretized using the consistent sample rate. 
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CHAPTER THREE 


EXPERIMENT 

An aeroelastic experiment was conducted in the Duke University Engineering wind 
tunnel facility. The goals of this test were to validate the analytical calculations of non- 
critical mode characteristics and to explicitly examine the aerodynamic mode divergence 
phenomenon. Additionally, because analyses show that the dynamic response of the 
system does not indicate divergence, a secondary goal of the testing was to evaluate 
different divergence onset prediction methodologies. To these ends, the simplest 
applicable model that could be devised was designed, fabricated and tested. 


Model Design 


The model design process first required that the desirable traits of the model be 
identified. Thus the non-dimensional parameter space was then examined to identify 
regions of parameters which would produce the best design. Physical parameter spaces 
were then examined to determine a configuration that could be built out ot realistic 
materials and tested in the facility available and with reasonable expectations ot 
instrumentation and data processing techniques. 


The desirable traits that were utilized in the model design range from the patently obvious 
to the sublime. It was desirable to have the model shaped like an airfoil. This is an 
important limitation- the shape restricts the strength, stiffness and inertia combinations 
which are achievable. It is also desirable that the model be constructed of machinable 
materials. This is a limitation, particularly in terms of an upper limit on the density and a 
lower limit on the strength and stiffness of available materials. It was also desired to 
have a wing structure that would not introduce additional modes into the experiment. 

This required that the airfoil be rigid in both the chord-wise and span-wise directions. 

The facility in which a model is tested places additional limitations on the design space. 
The model design must diverge at a dynamic pressure that the wind tunnel can reach. To 
ease mounting and eliminate tip aerodynamic phenomena, it was desirable to have a wing 
which would span the entire tunnel. To eliminate gravitational effects from the 
experiment, it was desirable to mount the wing between the floor and ceiling rather than 
spanning from one sidewall to the other. These two requirements fix the span. It is 
desirable to avoid tunnel blockage effects. The cross-sectional area of the tunnel, 
coupled with expected deflection angles of the model, sets an upper limit on the airfoil 
chord length if blockage is to be avoided. 
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The above requirements limit the design space before consideration of demonstrating the 
desired phenomenon are addressed. It is desirable to find a place in the design space 
which satisfies all of the above constraints and also demonstrates aerodynamic mode 
divergence. For clarity in identifying that the system has diverged while the dynamic 
mode persists, the designed model should have a high structural dynamic mode frequency 
at divergence. Low damping of the dynamic mode at the divergence condition is also 
desirable. The lower the damping of this mode, the more obvious the modal 
characteristics will appear. Further, it is desirable to demonstrate, through simple 
configuration changes to the model, different types of divergent behavior. 

The non-dimensional parametric variation database of modal characteristics at divergence 
was discussed previously. The structural dynamic parameters varied in the database are 
elastic axis location, e/b, radius of gyration, r u , and mass ratio, p. This database was 
utilized to identify regions in the parameter space where the divergence mechanism 
exhibited aerodynamic lag divergence. Additionally, it served to find regions where the 
damping and frequency characteristics were most desirable. 

Dimensional design spaces were also constructed. Although not presented here, the 
torsional spring stiffness, pitch inertia, mass, semichord and elastic axis location 
variations were examined. These variations allowed consideration of many of the issues 
associated with model construction and testing. They served to reduce the design space 
to a general description of the model. Fundamental characteristics of the configuration 
chosen were a NACA 0012 airfoil with a chord length of 8 inches, manufactured from 
aluminum as a thin-walled closed cell with spanwise stiffeners located near the elastic 
axis. Trailing edge segments made of different materials serve as the mechanism to 
determine and reconfigure the pitch inertia. 

Design-specific variations of physical variables were then examined. There were three 
physical quantities that were still adjustable within reasonable limits and still produce a 
model which could be manufactured: torsional spring stiffness, mass of the trailing edge 
segment, and distance from the trailing edge segment to the center of rotation. 

The torsional spring stiffness could be adjusted. As noted in the non-dimensional 
parametric variation, this produced no effect on the migration pattern of the eigenvalues. 
The torsional spring stiffness was used to control the magnitude of the frequency of the 
dynanic mode at the divergence condition. This is a measurement and data processing 
fidelity issue, not a phenomenon issue. The limitation on increasing the stiffness is the 
tunnel dynamic pressure capability. Increasing the stiffness increases the frequency, but 
also the divergence airspeed. 

The two remaining design freedoms were simultaneously varied; the results of these 
variations are shown in Figure 32 and Figure 33. These figures are very similar in shape 
and magnitude to the sample results presented for the non-dimensional parameter 
database at a fixed mass ratio, Figures 15 and 16. Figure 32 presents the ratio of the 
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dynamic mode frequency at the divergence condition to the air-off torsional frequency. 
The undulating surface is a function of elastic axis position and ratio of the trailing edge 
mass to the total mass of the system. The solid square shown at the ratio of masses value 
of 0.01 corresponds to the Plexiglass trailing edge configurations, configurations # 2 and 
# 3. The solid symbol at the ratio of masses value of .56 represents the Tungsten trailing 
edge configuration, configuration #1. From this figure, it is anticipated that the two 
designed inertia configurations will produce different dynamic mode migration. 


1 


-10.9 



Figure 32 Physical parameter variation results; ratio of frequency of dynamic mode at divergence to 
pitch mode natural frequency 

Figure 33 presents the damping of the dynamic mode at the divergence. Again, the 
surface is a function of elastic axis position and ratio of the trailing edge mass to the total 
mass of the system. 
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Figure 33 Physical parameter variation results; damping of dynamic mode at divergence 
condition 


Configuration Descriptions 

Three configurations of the model were designed for testing. The configurations differ in 
their torsional stiffness and inertial properties. These properties influence the natural or 
zero-airspeed structural dynamic frequency. In addition, analyses indicated that the 
subcritical migration pattern of the dynamic mode eigenvalue is strongly influenced by 
the inertial properties. 

Reconfiguring the model’s pitch inertia was accomplished by changing the airfoil trailing 
edge. This also changes the non-dimensional mass ratio. Configuration #1 employed the 
Tungsten trailing edge component, while the second and third configurations employed 
the Plexiglass trailing edge component. These materials were chosen to provide a large 
difference in the torsional inertias, and thus an observable difference in the non-critical 
mode behavior. Examining Figure 32 and Figure 33, there are solid symbols shown on 
the surfaces at the ratio of masses value of .56; these are the expected frequency ratio and 
damping for the Tungsten trailing edge configuration, configuration #1 . The solid 
squares shown in the figures at the ratio of masses value of 0.01 corresponds to the 
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Plexiglass trailing edge configurations, configurations # 2 and # 3. The two inertia 
designs lie on different tiers of the frequency ratio surface. At divergence, the frequency 
of the tungsten configuration is expected to drop to 29% from of its air off frequency. 
The frequency of the Plexiglass trailing edge configuration is expected to drop to 53% of 
the air off frequency when the system destabilizes. Additionally, the dampings are 
different. The analytical frequency ratios and damping at divergence for the three 
configurations are given in Table 1 1 . From this data, it is anticipated that models with 
the two designed inertias will produce discernibly different dynamic mode migrations. 


Config 

mass T railingEdge 

Wa 

g>d 

®d/ 


# 

mass Total 

(air off) 


/®a 


1 

0.56 

21.2 

6.2 

0.29 

0.51 

2 

0.01 

49.5 

26.4 

0.53 

0.75 

3 

0.01 

87.3 

46.4 

0.53 

0.75 


Table 1 1 Analytical frequency ratios and damping at divergence 

Reconfiguring the model stiffness was accomplished by changing the torsional spring and 
thus the stiffness. The first two configurations used a 1-inch diameter torsional spring 
with an advertised stiffness of 0.94 lb r in/degree. The third configuration employed a %- 
inch spring tine with an advertised stiffness of 3.18 lbf-in/deg. The stiffnesses were 
measured and will be discussed in the experimental results section. 

Table 12 provides a summary of the configurations. The dimensional quantities and non- 
dimensional parameters for each were previously listed in Table 5. Model configuration 
#2 serves as a comparison configuration to each of the others and will be discussed in 


much more detai 

than the other two. 


Trailing Edge Segment 

Torsional Spring 

Config # 

Material 

masSfE / 

/ mass Total 

Diameter 

Stiffness 

(lbf-in/deg) 

1 

Tungsten 

0.56 1 

1” 

0.90 

2 

Plexiglass 

0.01 

1” 

0.90 

3 

Plexiglass 

0.01 


2.78 


Table 12 Description of wind tunnel model configurations 
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Hardware 


Wind Tunnel Model 

The divergence assessment testbed ( dat ) wind tunnel model consists of a typical section 
airfoil with a flexible mount system providing a single degree of freedom structural 
dynamic mode. The only structural dynamic mode of this model is torsional rotation, or 
angle of attack. 


Airfoil Description 

The airfoil section is a NACA 0012 with an 8-inch chord and a span of 21 inches. This 
spans the entire test section from the floor to ceiling, as shown in Figure 34. The airfoil 
is an aluminum shell, 1/32 inch thick. To ease fabrication and instrumentation it was 
made in two sections that join at approximately the mid-span. The internal structure has 
two spar webs running the entire span to provide bending rigidity and the designed 
inertial properties. Each of the two span sections consists of internal spars and airfoil 
which were cut as a single entity from a solid block of aluminum using a wire electro- 
deposit-machine (EDM). The last 1.125 inches of the airfoil were fabricated separately to 
provide test configurations with different inertial properties. To effect a large change in 
inertia, trailing edge segments were fabricated from Plexiglass and from Tungsten. These 
trailing edges could be easily changed during the test. 


67 



Figure 34 Researcher installing wind tunnel model; airfoil shown with Plexiglass trailing edge 
segment 

Mount System 

The mounting system for the dat model has a ceiling mechanism and a floor mechanism. 
Both portions of the mount system are required for mounting the model and holding it in 
place. 

The ceiling mechanism, shown in Figure 35, serves three functions in addition to holding 
the model in place. The torsional spring is contained in the ceiling mechanism. 
Sometimes called a barrel spring or a Bendix flexure, the spring provides the stiffness 
associated with the structural dynamic torsion mode. The ceiling mount also contains a 
turntable which allows the rigid angle of attack to be set and changed. Mounted between 
the turntable and the torsional spring is a balance which measures torsional strain. 



Figure 35 Ceiling mechanism of mount system 
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This mount system was designed and built with couples to connect the torsional spring to 
the balance and to the airfoil. Three sets of couplers were fabricated so that springs of 
various sizes could be easily substituted. This system flexibility allows configuration 
changes in stiffness during the test. 

The floor mount system serves two functions in addition to holding the model in place. 
The floor mount mechanism is shown from below, looking up at the airfoil and ceiling 
mount system. Figure 36. 



Figure 36 Floor mount mechanism, viewed from below 

A shaft extends out of the wing through the tunnel floor and is fixed to an angular 
displacement transducer. This allows the total angle of attack of the wing to be measured 
including the rigid angle of attack set using the ceiling mount and the elastic increment. 
Additionally, the sensor wing and pressure reference tubes pass out of the tunnel through 
this mount system. 


Torsional stiffness testing of springs 

The torsional springs were tested to determine their stiffness constants, as well as 
evaluate the range of operation and linearity. The springs were inserted into a test fixture 
which measured the torsional moment and the deflection angle. Weight was applied to 
the test fixture in a manner which caused one end of the spring to rotate through a 
deflection angle. The applied weight was increased and data acquired at each weighting 
until the deflection angle stayed constant. The data set was fit with a linear equation that 
minimized the error in the least squares sense. The torsional spring stiffness constant is 
the slope of torsional moment versus deflection angle. 

The data for the l-inch diameter spring is shown in Figure 37, along with a linear fit to 
the data and the manufacturer’s specifications. Several data sets were acquired for each 
spring. Each data set was curve fit with and without constraining the y-intercept. Curve 
fits to the data sets produced spring constants between .87 and 1.08 lbrin/degrees. After 
consideration of pluck test results and inertia measurements, the value to be used in 
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comparison analyses was set as 0.90 lbrin/degree. The manufacturer’s specification for 
this spring, 0.94 lbrin/degree, falls within the scatter of the measurements and curve fits. 



Angle (degrees) 

Figure 37 Stiffness data for 1 inch diameter torsional spring 


The data for the % -inch diameter spring is shown in Figure 38, along with a linear fit to 
the data and the manufacturer’s specifications. Curve fits to the data sets produced spring 
constants between 2.6 and 2.8 lbrin/degrees. The nominal value for performing 
comparison analyses was set as 2.78 lbrin/degree. The manufacturer’s specification for 
this spring, 3.18 lbrin/degree, is significantly higher than the measured value. This 
difference invalidated preliminary estimates of the divergence dynamic pressure. 
However, the qualitative behavior of the eigenvalues were shown in the analysis to be 
invariant with respect to torsional stiffness. 

Based on the measured data, both springs appear to behave in a linear manner for angular 
displacements below 18 degrees. 
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Figure 38 Stiffness data for 3 A inch diameter torsional spring 


Model Instrumentation 

The wind tunnel model was instrumented with the following sensors. A Trans-tek 
angular displacement transducer was mounted to the bottom shaft of the airfoil. The 
measurement was made such that the total angle of attack was measured, including the 
prescribed rigid contribution. A torsional moment sensor connected the angle ot attack 
turntable to the torsional spring coupler. This sensor served the role traditionally played 
by a balance. A piezoresistive accelerometer was mounted inside the airfoil towards the 
leading edge. This type of accelerometer is dc -coupled, providing static values, as well 
as dynamic frequency responses. Unsteady pressure transducers were mounted at the 
approximate midspan of the airfoil section. One transducer was on the leading edge. 
Eight others comprise pairs of top- and bottom- mounted pairs. These eight transducers 
were mounted at the 5%, 15%, 30% and 52% chord locations, respectively. 

Facility Description 

This test was conducted in the wind tunnel facility at Duke University. Two gust vanes 
were installed vertically in the tunnel, ahead of the model. Each gust vane system 
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consists of a fixed airfoil and a slotted cylinder. The slotted cylinders were placed just aft 
of the gust vane airfoils. Each cylinder was mounted to a shaft. The shaft was held by 
rotating bearings in the ceiling and floor. Beneath the wind tunnel floor, the shafts were 
connected by a gear and belt system to a motor. Via this motor, the cylinders were 
rotated. Voltage supplied to the motor governs the motor speed and thus the rotational 
frequency of the slotted cylinder. Through this mechanism, aerodynamic forcing of the 
test article was achieved. 


Data Acquisition System 

The data was acquired using a portable, self-contained system. The signals were input to 
signal conditioners, the type of which depends on the sensor. The signal conditioners 
sent voltages to a National Instruments DAQ 700 card which was inserted into the game 
card slot of a laptop computer. 

The DAQCard-700 is a Type II PCMCIA card with 12-bit analog to digital conversion. 1 
It can acquire data from 1 6 single-ended analog input channels with a maximum sample 
rate of the card is 100 kilosamples per second. The card has DC input coupling, enabling 
static measurements. A FIFO buffer contains the data during multiple analog to digital 
conversions to prevent data loss. An onboard counter/timer generates the sample interval 
clock. For this study, multichannel acquisition with continuous scanning was the mode 
of operation, taking one reading per sample interval, always in the same channel order, 
starting with the highest channel number. 

The data acquisition system was driven by virtual instruments designed in Labview, a 
National Instruments software interface to their data acquisition products. 


Experimental Data Acquired 

Wind tunnel test data was acquired to investigate experimentally system behavior and 
validate analytical results. Specifically, data was acquired: 1) to find the divergence 
dynamic pressure; 2) to examine the modal characteristics of non-critical modes, both 
subcritically and at the divergence condition; 3) to examine the eigenvector behavior. A 
secondary goal of acquiring and analyzing this data was to evaluate standard 
experimental divergence onset prediction methods. Addressing these goals required that 
several different experimental methods be employed for acquiring data. 


1 DAQCard-700 User's Manual: Multifunction I/O Board for the PCMCIA BUS 
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The first type of data acquired was system response to turbulence occurring naturally in 
the tunnel. Initial divergence testing was performed by increasing the dynamic pressure 
or velocity of the airstream. Stability or instability of the system was observed; data was 
acquired. The majority of the turbulence-excited data was acquired at fixed wind tunnel 
velocity. These data are referred to in this text as velocity stabilized points. 

Forced response data was also acquired. One method employed to excite the system was 
applying a pluck or a rap to the mount system, as shown in Figure 39. A portion ot the 
mount system accessible from outside the tunnel was connected to the moving portion of 
the torsional spring. Applying a force to the bottom coupler was equivalent to rapping 
the airfoil. Data was acquired as the model was plucked or rapped; the model response 
was recorded as the motion decayed, or grew in the case of unstable systems. 



Figure 39 Administering a pluck to the model 


The second method of forcing the system was to employ the gust vanes, which were 
described previously. While the pluck test method utilizes the structure to apply the 
forcing function, the gust method acts through aerodynamic forces. Two types of gust 
excitations were used to acquire data. The first type of excitation used was a sweep of 
the frequency range of interest, starting at the high frequency end and progressing to low 
frequency. The second type of excitation provided by the gust vanes was rotation of the 
cylinder at a fixed frequency, dwelling at a single frequency. The latter method is 
referred to in this text as the sine dwell method of excitation. 

Data Processing 

The measured data was utilized to glean insights into the system behavior. Examining 
different system characteristics required that the data be processed in several ways. 

Identifying the divergence condition was accomplished by observing the system stability 
as the dynamic pressure was increased. Additionally, the divergence onset was predicted 
using several methods which are classically employed in experiments. Because it is a 
contention of this work that the dynamic mode persists at a nonzero frequency as the 
instability is reached, it is interesting to contrast the results obtained from static and 
dynamic prediction techniques. The static methods used were load monitoring, angle of 
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attack ratio monitoring, and the Southwell method. The dynamic methods employed 
were frequency tracking and inverse amplitude tracking. 

Identifying modal characteristics was accomplished by scrutinizing the time history data. 
Frequency domain analyses and the logarithmic decrement technique of determining 
damping were utilized. The frequency domain analyses include analyzing the Fourier 
transformations, the power spectral density functions and approximated transfer 
functions. 


Experimental Results 
Configuration #2 

The experimental data for wind tunnel model configuration # 2 is presented first. 
Configuration # 2 has the l-inch diameter torsional spring and has the Plexiglass trailing 
edge. This configuration has lower torsional inertia than the Tungsten configuration, 
model configuration #1, but has the same torsional stiffness. This configuration has 
lower torsional stiffness than the % inch spring configuration, model configuration # 3, 
but has the same torsional inertia. Thus this configuration serves as a comparison 
configuration to each of the others. Results will be presented and discussed in the 
following order: determination of the divergence condition, subcritical techniques for 
predicting divergence onset, system behavior at divergence, and subcritical modal 
characteristics. 


Divergence dynamic pressure 

Divergence testing was conducted by setting the zero airspeed angle of attack, Oto, which 
is referred to as the rigid angle of attack, as close to zero as possible. The divergence 
dynamic pressure was determined by increasing the velocity and recording data as the 
system became unstable. A time history showing the divergence of this configuration is 
shown in Figure 40. The dynamic pressure was being slowly increased until the angle of 
attack increased dramatically and suddenly. This was declared as the divergence 
dynamic pressure, 5. 1 psf (244 N/m 2 ). The time history shows that the model oscillates 
about a new angle of attack position, which is not at the hard stop of the spring. It is 
speculated that the airfoil has reached an angle of attack where flow has separated and 
stall has occurred. This issue and the ensuing behavior of the system are further 
discussed at the end of this section on experimental results. 


74 




Figure 40 Divergence of wind tunnel model configuration #2 


Divergence prediction using subcritical data 

Five methods of experimentally predicting divergence onset were utilized in this work. 
Three methods which examine the static properties and two dynamic response methods of 
analyzing data were applied. 

Static load monitoring is a fundamental method of predicting divergence. This method 
depends on the monitored load becoming large as divergence is approached. The slope 
of the moment versus dynamic pressure curve is the key parameter. This slope changes 
dramatically in the neighborhood of divergence. In this experiment, the torsional spring 
moment is the monitored load. In applying this method, data sets are acquired at several 
rigid angles of attack. For each rigid angle of attack, data was recorded at regular 
intervals of dynamic pressure. Data for each angle is treated as a data set. The data is 
shown in Figure 42. The legend indicates the rigid angle of attack for each of the traces 
in the figure. 
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Two data sets were acquired for rigid angles of attack very near zero degrees as shown in 
Figure 41 . For these data sets, the torsional moment is very small until just before the 
divergence condition is reached; the load is seen to increase dramatically at this point. 
The steep slope indicates that divergence is imminent. 

Data is also presented for non-zero rigid angles of attack. The terminations points of 
these curves show the last dynamic pressure before the system became unstable for each 
rigid angle of attack. A method for estimating these destabilizing dynamic pressures will 
be presented in the discussion section. The increase in load with changing dynamic 
pressure is more gradual for larger a () . The structural moment is directly proportional to 
the elastic increment in the angle of attack. Examination of the static equilibriun 
equations for this system shows that the elastic increment is an amplification of ot«. 



a 


-0.0305 ! 

0.0707 

1.1405 

2.1057 

3.0114 

4.0286 

4.8911 

6.0194 


Figure 41 Strain monitoring for predicting divergence onset 


The second method of predicting the onset of divergence is examining the angle of attack 
as the dynamic pressure is increased. Divergence is classically defined by the angle of 
attack becoming infinitely large, or the inverse of the angle becoming zero. Figure 42 
shows the inverse angles of attack data. The data has been normalized by the rigid angle; 
this normalization collapses the data to a single line. 
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Figure 42 Predicting divergence using angle of attack ratio 

A linear fit through the experimental data produces a divergence dynamic pressure of 5 
psf. The theoretical curve, which employs the measured value of torsional stiffness 
produces a divergence dynamic pressure of 4.6 psf. The experimental prediction using 
this method is 9% higher than the theoretical result. The disagreement between analysis 
and experiment is resolved by utilizing a measured value of lift coefficient instead of the 
ideal of 2k. Using Cux = 5.7 1 /radians produces a divergence dynamic pressure of 5.08 
psf. 

The two static methods used previously are combined in the Southwell method. To apply 
the method, the static load is measured at fixed dynamic pressures for various rigid 
angles of attack. The data at each dynamic pressure constitutes a single data set. A linear 
fit is made to each data set, plotting static moment versus angle of attack. The data is 
shown in this manner in Figure 43. The slope of each line is denoted, X. Divergence 
occurs at the dynamic pressure which makes the slope of these data infinite. 
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Figure 43 Intermediate Southwell plot 


The data sets are combined into a single plot. The slope, X, is plotted versus itself, 
normalized by the dynamic pressure. The slope of this line predicts the divergence 
dynamic pressure. Figure 44 shows the data for this configuration, with a slope of 5.5 
psf. 



Figure 44 Southwell method for predicting divergence onset 


During a divergence onset test, the Southwell method is applied as the dynamic pressure 
data is acquired. Figure 45 shows the prediction of the Southwell method as additional 


78 




data is considered in the linear fit. Using only data below 2. 1 psf, the divergence 
prediction of 5.2 is fairly close to the final divergence prediction which utilized data up to 
5.2 psf. 
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Figure 45 Southwell method results using increasing amounts of data 

Dynamic methods were also applied to predict divergence. Divergence is classically 
considered to occur as the torsional mode frequency drops to zero and then statically 
destabilizes. One classical method of predicting the onset of divergence is to monitor the 
torsional mode frequency migration, anticipating that it will go to zero prior to 
divergence. Figure 46 shows the system frequency extracted from subcritical data as the 
dynamic pressure is varied. The Fourier transformation of the angle of attack response 
was computed for data generated when frequency sweeps were input to the gust vanes. 

At divergence, the frequency is still greater than 3 Hz for this configuration. 
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Figure 46 Frequency tracking to predict divergence 


A second method that relies on the same data and analytical techniques as the frequency 
tracking method. The amplitude of the system response is anticipated to become large as 
divergence is encountered. Rather than utilizing the static response as was done 
previously in the load monitoring method, the amplitude of the modal response is 
tracked. Figure 47 shows this data. The amplitude is actually seen to decrease until the 
last data point before divergence. 



Figure 47 Inverse amplitude of the power spectral density of the angle of attack response 
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Extracting modal characteristics 


One of the primary concerns of this study is tracking the modal characteristics of the 
aeroelastic system as the dynamic pressure increases towards the divergence condition. 
Identification of the torsional mode frequency and damping was accomplished by 
employing different excitation, measurement and data processing techniques. 

The structural damping of the system without any aerodynamic forces acting can he 
determined by considering the decay of the response due to a rap or a pluck. A pluck was 
applied to the model, as shown in Figure 39 and the logarithmic decrement method was 
used to analyze the data. 

The logarithmic decrement is a well-documented method for calculating the damping of a 
system 2 . Applying the method is dependent upon having a time history of a decaying 
response where several cycles of the decaying motion are evident. Air off pluck test data 
is shown in Figure 48. The time history shows the angle of attack displacement as a 
pluck is administered to the model and the motion decays. The circular symbols indicate 
the values which were used in the damping calculation. 



Time, sec 


Figure 48 Extracting logarithmic decrement data from angle of attack time history, air off data for 
configuration # 2 

The logarithmic decrement, 5, is calculated using Equation 20 which then produces a 
damping value from Equation 2 1 . The amplitudes of the 0 th and N th cycles are denoted 

2 Clough. Ray W., and Joseph Penzien, Dynamics of Structures, 
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by a in this equation. The real part of an eigenvalue can be calculated from the damping 
and the frequency, as shown in Equation 22. 




log 


«<) 


N 


cycles 


Equation 20 




/vW 2 + s 2 


Equation 21 


V “ ~ Equation 22 

l-sin 2 (f) 


Data was acquired and analyzed in this manner for subcritical dynamic pressures. Figure 
49 presents the results of numerous applications of the logarithmic decrement method. 
The damping and frequency information were converted to real and imaginary 
components representative of a single mode. The real part is plotted as a function of the 
dynamic pressure. Results from individual time histories are shown as small open 
circles. The larger solid circles are the average results at each dynamic pressure. While 
the damping is well defined and repeatable at low dynamic pressures, the uncertainty at 
the higher dynamic pressures becomes quite significant. 
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Figure 49 Logarithmic decrement results; large symbols are average values for each 
dynamic pressure 


Frequency domain analysis is an important tool in identifying the modal characteristics of 
a system. By performing a Fourier transformation on a set of data, dominant frequencies 
of the system become evident. Damping can also be gleaned from data examined in this 
manner. The frequency of the system is determined by the maximum amplitude. 
Extracting the damping requires that one find the amplitude of the peak and then the 

frequencies where the amplitude is reduced by a factor of V^, these frequencies are 
referred to as the half power points. Figure 50 illustrates the quantities required to 
perform these calculations. This method was applied when data across the frequency 
range of concern could be obtained in a single time history, as in the case of the 
frequency sweep excitations. 
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Figure 50 Frequency domain representations for extracting modal parameter data 

Figure 5 1 shows the Fourier transformation of the angle of attack data acquired at 1 .52 
psf for configuration # 2. This data was produced by exciting the system with a 
frequency sweep signal to the gust vanes. From this plot, the maximum amplitude was 
found to be 766 at a frequency of 7.3 (45.6 rads/sec). The half power points were found 
to be 6.75 Hz and 7.6 Hz. The difference between these frequencies, denoted Ato, is 0.85 
Hz. The damping is calculated using Equation 23, resulting in a value of C, = 0.059, or an 
eigenvalue with real part -2.67 . 


A oj 
2cq 


Equation 23 


This half power method was also utilized in real-time. A sine dwell excitation was sent 
to the gust vanes. The frequency of the excitation was tuned until the center frequency 
was determined. The amplitude was recorded and the half amplitude was determined. 
The frequency was then retuned to determine the half power points. In addition to 
recording these important values, time histories were recorded as the model was excited 
at each of these three frequencies. When these time histories were processed, they were 
found to be slightly inconsistent with the real-time determination of the half power 
points. The recorded and processed data was utilized in the final computations pertinent 
to the real-time half power method. 
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A variation of the Fourier transformation method uses the power spectral density 
function, PSD. The PSD is the normalized square of the Fourier transformation. The 
normalization allows frequency data extracted from different time length and sample rate 
data to be compared. Because the PSD is essentially the squared frequency content, the 
half power points are assessed at the peak amplitude/2. The right side of Figure 50 
illustrates the required quantities to apply this technique. It is identical in calculation to 
the Fourier transform method. 

This method was applied for velocities in the vicinity of divergence. Very near the 
dynamic pressure at which the system destabilized, sine dwell excitations had to be used 
to impart enough energy to the dynamic mode. Analysis of each time history produced a 
single point on a frequency domain plot. 



Frequency (Hz) 


Figure 51 Amplitude of Fourier transformation of angle of attack; Configuration # 2, 
dynamic pressure 1.52 psf 


The transfer function is another frequency domain interpretation which can be utilized to 
extract modal information. A transfer function requires an input, or excitation, signal in 
addition to the response signal. The actual excitation, the pressure field induced by the 
gust vanes on the wing is not available. The dynamic pressure of the tunnel, however, 
can be utilized in some sense as the excitation. Ewins 3 suggests plotting the real part of 
the transfer function as a function of frequency. This is done using sine dwell data to 
produce each point. 


’ Ewins, D. J., Modal Testing: Theory and Practice 




The objective of the logarithmic decrement analysis and the frequency domain analyses 
were to determine the subcritical modal characteristics of the system as the dynamic 
pressure was increased towards divergence. The results of these techniques, applied to 
configuration #2, are summarized in Figures 52, 53 and 54. The damping and frequency 
information have been converted to the real and imaginary parts of an eigenvalue 
assuming that they were representative of a single mode. 

Figure 52 shows the real part as a function of dynamic pressure for the different methods. 
This is indicative of the damping of the system. Seven sets of data are shown in the 
figure. All data sets indicate that the damping of the system is mainly due to the 
aerodynamics. The structural damping, indicated near zero airspeed is very small by 
comparison. All methods show the modal damping increasing as the dynamic pressure 
increases until divergence is reached. The individual curves are discussed below. 
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Figure 52 Real part of measured modal data as functions of dynamic pressure, 
Configuration # 2 

The frequency sweep data was processed using the Fourier transform method. Because 
the frequency sweeps did not impart enough modal energy to the model near divergence, 
the results are only shown for dynamic pressures below 4. 1 psf. These subcritical results 
show a monotonic increase in the magnitude of the real part as dynamic pressure 
increases. The data labeled as “maximum damping” corresponds to using the outside 
edges of the modal peak to determine the half power frequencies. The data labeled “best 
guess” corresponds to faring a curve through the scatter of the frequency data. The line 
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faring is approximately equivalent to processing a shorter time segment, which increases 
the Fourier transform frequency spacing. The “best guess” produced lowest values of 
damping across the dynamic pressure range. These values, which were more subjectively 
determined than the maximum damping values, produced more variability in the results. 
The “maximum damping” interpretation of the frequency sweep data lies nearly on top of 
the logarithmic decrement results. 

The logarithmic decrement results, shown in the figure by open squares, are the average 
values at each dynamic pressure, which were shown in Figure 49. The logarithmic 
decrement results are the only results available for the air off condition. The results show 
very low structural damping; for this configuration it was measured as £ = 0.0053. 
Translated into the real part of an eigenvalue whose frequency is 49.5 radians/second, 
this corresponds to a real part of -0.26. A very orderly march to higher damping is 
charted by this data, which is available only at dynamic pressures below 3.6 psf. 

The largest damping values were obtained by the real-time half power method. This data 
was taken up to 4.6 psf and provides a very smooth trend in the damping values. 

In the vicinity of divergence, the sine dwell data was analyzed using the frequency 
domain methods. These analyses provide interesting trends. The Fourier transform 
results, shown by open triangles connected by a solid line, begin near where the 
frequency sweep data ended. The data continue to decrease until a dynamic pressure of 5 
psf is exceeded. The values decrease as dynamic pressure is further increased. The 
power spectral density method and the transfer function method yield similar trends. 

The frequency of the dynamic mode is shown in Figure 53. All methods show the same 
trend. The structural dynamic frequency, determined from time histories is 49.5 
radians/second, (7.9 Hz). The frequency decreases as the dynamic pressure and thus the 
aerodynamic coupling increase. All results lie on top of each other except for the 
analyses of the sine dwell data. These data were acquired near divergence and have some 
scatter in the results. Regardless of the method, however, the frequency of the dynamic 
mode is determined to be not lower than 20 radians/second (3.2 Hz) at the divergence 
condition. The dynamic mode is clearly shown to persist at a non-zero frequency at the 
divergence condition. 
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Figure 53 Imaginary part of measured modal data as functions of dynamic pressure. 
Configuration # 2 

The subcritical modal data is also presented in a root locus format. Figure 54. The sine 
dwell results have been removed from this chart for purposes of clarity. The four traces 
shown, which are not individually identified, correspond to the pluck test results analyzed 
using the logarithmic decrement technique, the real-time half power results and the two 
interpretations of the frequency sweep data analyzed with the Fourier transform. All of 
the traces begin in the neighborhood of the air-off value at -0.26 + j49.5 and migrate 
generally downward and to the left as dynamic pressure increases. 
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Configuration #1 


The experimental data for wind tunnel model configuration # 1 is now presented. 
Configuration # 1 has the 1-inch diameter torsional spring and has the Tungsten trailing 
edge. The Tungsten trailing edge produces a higher torsional inertia than for the other 
configurations. 

Divergence testing was conducted exactly as for configuration #2. The zero airspeed 
angle of attack was set as close to zero as possible. The divergence dynamic pressure 
was increased slowly until the angle of attack increased dramatically and suddenly, as 
shown in the time history of Figure 55. This was declared as the divergence dynamic 
pressure, 5.14 psf (246 N/m 2 ), which is nearly identical to the divergence condition for 
configuration #2. Based on analysis, this was the anticipated result. 



Time (seconds) 


Figure 55 Divergence of wind tunnel model configuration #1 

The subcritical methods of predicting divergence onset were each applied for this 
configuration. The data are similar to that presented for configuration #2 and are omitted 
here. The Southwell method was shown to be the most reliable prediction technique in 
analyzing the data for configuration #2. For configuration #1, a divergence condition of 
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5. 1 psf was determined. This prediction required that data be known up to 3. 1 3 psf, or 
61% of the predicted divergence dynamic pressure. Using data at a higher dynamic 
pressure was found to detract from the accuracy of the prediction, as nonlinear effects 
began to influence the data. 

The subcritical modal characteristics of this configuration were determined, as for the 
previous configuration, as the dynamic pressure was increased towards divergence. The 
results of these techniques, applied to configuration #1, are summarized in Figures 56, 57 
and 58. The damping and frequency information have been converted to the real and 
imaginary assuming that they were representative of a single mode. Only two methods 
were used to acquire and analyze the data for this configuration: the logarithmic 
decrement technique and the half power method applied to frequency sweeps. 

Figure 56 shows the real part as a function of dynamic pressure. As with configuration 
#2, all data sets indicate that the damping of the system is mainly due to the 
aerodynamics. The logarithmic decrement results are shown in. the figures by open 
squares. The air off results show very low structural damping; for this configuration it 
was measured as ^ = 0.0046, a very low value by comparison. The half power method 
results are shown in the figures by solid diamonds. Both methods show the modal 
damping increasing as the dynamic pressure increases for subcritical conditions below 4 
psf. The measured damping tends towards zero near the divergence dynamic pressure. 
This issue is further addressed in the discussion section of this paper. 
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[ □ Pluck Test Data, Log Decrement Method 
♦ Frequency Sweep Data, Half Power Method 
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Figure 56 Real part of measured eigenvalues. Configuration # 1 

The frequency of the dynamic mode is shown in Figure 57. The structural dynamic 
frequency of the air off pitch mode was determined as 21.2 radians/second, (3.4 Hz). The 
logarithmic decrement and half power method show the same trends for increasing 
airspeed . The frequency decreases as the dynamic pressure and thus the aerodynamic 
coupling increase. The frequency of the dynamic mode is determined to be not lower 
than 7 radians/second (1.25 Hz) at the divergence condition. The dynamic mode is 
clearly shown to persist at a non-zero frequency at the divergence condition. 
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Figure 57 Imaginary part of measured eigenvalues, Configuration # 1 

The subcritical modal data is also presented in a root locus format. Figure 58. The traces 
shown correspond to the pluck test results analyzed using the logarithmic decrement 
technique and the frequency sweep data analyzed with the Fourier transform. All of the 
traces begin in the neighborhood of the air-off value at -0. 1 + j21.2 and migrate generally 
downward and to the left as dynamic pressure increases. Near divergence, they resemble 
a left-to-right scribbling due to the uncertainty in determining the damping value. 
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Figure 58 Root locus using experimental data. Configuration # 1 


Configuration #3 

The experimental data for wind tunnel model configuration # 3 is presented next. 
Configuration # 3 has the % -inch diameter torsional spring and has the Plexiglass trailing 
edge. The inertial characteristics are identical to configuration #2. The torsional stiffness 
has been changed. As previously mentioned, changing the torsional stiffness is not 
anticipated to affect the dynamic characteristics of the system; in the anlaysis, the 
eigenvalue migration pattern remains identical. Although the nondimensional root locus 
remains the same, the dimensional natural frequency of the pitch mode increases as does 
the dynamic pressure which produces divergence. 

Divergence testing was conducted exactly as for configuration #2. The zero airspeed 
angle of attack was set as close to zero as possible. The divergence dynamic pressure 
was increased slowly until the angle of attack increased dramatically and suddenly, as 
shown by the angle of attack time history. Figure 59. This was declared as the 
divergence dynamic pressure, 15.2 psf (730 N/m~). 
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Figure 59 Divergence of wind tunnel model configuration #3 

The subcritical methods of predicting divergence onset were each applied for this 
configuration. For this configuration, the Southwell method predicted a divergence 
condition of 15.8 psf. Although the prediction is fairly close to the divergence dynamic 
pressure observed directly, the nonlinear effects influenced this data more significantly 
than the previous configurations. The linear approximation to the data, utilized in the 
Southwell method, did not fit the data well. As the range of dynamic pressures included 
in the data set analyzed increased, the predicted divergence dynamic pressure also 
increased. 

The subcritical modal characteristics of this configuration were determined, as for the 
previous configurations, as the dynamic pressure was increased towards divergence. The 
results of these techniques, applied to configuration #3, are summarized in Figures 60, 6 1 
and 62. The damping and frequency information have been converted to the real and 
imaginary assuming that they were representative of a single mode. The damping and 
frequency traces presented are computed using the methods described pertaining to 
configuration #2. 

Figure 60 shows the real part as a function of dynamic pressure. As with configuration 
#2, all data sets indicate that the damping of the system is mainly due to the 
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aerodynamics. The logarithmic decrement results are shown in the figures by open 
squares. The air off results show very low structural damping; for this configuration it 
was measured as £ = 0.0035. Results from the other various methods are shown by the 
symbols indicated in the legends. All methods show the modal damping increasing as the 
dynamic pressure increases. Additionally, the data show that the dynamic mode is stable 
throughout the range of dynamic pressure, including at the divergence condition. 
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Figure 60 Real part of measured eigenvalues, Configuration # 3 

The measured frequency of the dynamic mode is shown in Figure 61. The structural 
dynamic frequency of the air off pitch mode was determined as 87.15 radians/second, 
(13.9 Hz). All methods show the same trends for increasing airspeed. The frequency 
decreases as the dynamic pressure and thus the aerodynamic coupling increase. At 
divergence, the frequency of the dynamic mode is determined to be not lower than 36 
radians/second (5.8 Hz) at the divergence condition. The dynamic mode is clearly shown 
to persist at a non-zero frequency at the divergence condition. 
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61 Imaginary part of measured eigenvalues. Configuration # 3 
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Figure 62 Root locus as a function of dynamic pressure, Configuration # 3 
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The subcritical modal data is also presented in a root locus format. Figure 62. All of the 
traces begin in the neighborhood of the air-off pitch mode eigenvalue, -0.3 + j87. 1 , and 
migrate generally downward and to the left as dynamic pressure increases. The trend 
shown in this figure is for the dynamic mode damping to increase dramatically in the 
neighborhood of the divergence dynamic pressure. 

Hard Limit Instabilities 

Time histories of each configuration at their respective divergence dynamic pressures 
have been presented, Figures 40, 55, and 59. Each figure shows that divergence resulted 
in the model sitting at an angle of attack where the airfoil has stalled. Data was acquired 
for velocities beyond these divergence conditions until further destabilization resulted in 
the model hitting the hard stops of the torsional spring. An example of this is shown for 
configuration #2 in Figure 63; the angle of attack time histories is presented. Beginning 
at 15.8 seconds into the time history, the character of the motion changes dramatically. 
Destabilization, characterized by the onset of this dramatic motion, occurred at a dynamic 
pressure of 5.59 psf (268 N/m 2 ). 



Figure 63 Hard limit instability encountered as dynamic pressure is increased, 
(Configuration #2) 

A second example is shown in Figure 64; this data was acquired while the tunnel 
condition is held constant at 5.53 psf. This data shown in both time histories indicates 
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that the system destabilizes in a static sense, while dynamic motion persists. 
Examination of the data using an expanded scale. Figure 65, shows the dynamic 
oscillations of the torsional motion at 4.9 Hz as the system becomes statically unstable. 



Figure 64 Divergence encountered while at constant dynamic pressure 



Time (seconds) 

Figure 65 Divergence encountered while at constant dynamic pressure, expanded scale 
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The time histories presented above appear to have different character, despite the same 
configuration being at approximately the same condition. The system behavior has to be 
examined in the light of nonlinearities and regions of aerodynamic and structural loading 
and unloading. Figure 66 addresses this issue. Shown in the figure are regions where 
structural and aerodynamic nonlinearities affect the system. The torsional springs have 
been shown to have limits on the linear behavior at approximately 18°. This deflection 
limit is relative to the neutral position of the spring. Because the rigid angle of attack was 
set by rotating both ends of the barrel spring, the neutral position is equivalent to the rigid 
angle of attack, oto. All structural loading or unloading is relative to the neutral position 
of the spring. The system also appears to be subject to an aerodynamic nonlinearity as 
the flow separates and stall is encountered. Stall initiates at an angle substantially below 
the structural limit; effects were observed as low as 8°. The stalled airfoil stabilizes at 
angles as high as 1 1 Vi°. Aerodynamic loading or unloading is relative to the symmetric 
position of the wing, 0°. These regions determine and explain the behavior of the system 
observed in the time histories. 
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Figure 66 Regions of behavior for the airfoil 

Consider the time history shown in Figure 55. Prior to the destabilization, shown to 
begin at 15.7 seconds, the linear system has already diverged. Because the airfoil is 
sitting at approximately 1 1 V 2 0 , aerodynamic stall has effectively decreased the lift curve 
slope; in a simplistic static sense, this can be thought of as a decrease in (negative) 
aerodynamic stiffness. A decrease in negative aerodynamic stiffness is raises the 
dynamic pressure at which the system destabilizes. Thus the nonlinear system appears 
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stable at a dynamic pressure which exceeds the linear system divergence condition. At 
15.7 seconds, the system possesses enough aerodynamic energy such that the nonlinear 
system becomes unstable. Because the airfoil is stalled only as it oscillates in one 
direction (in this case, toward larger angles), the motion is forced towards zero. As the 
system destabilizes, the airfoil reenters the linear aerodynamic range. The airfoil has 
enough momentum to approach zero, which is the aerodynamical ly and structurally 
unloaded position, since the airfoil is at zero degrees rigid angle of attack. The system is 
nowat an unstable equilibrium point. Pursuant motion may be to either the positive or the 
negative side. In the case of Figure 55, the linear range system diverges in the same 
direction that it just came from. In the case of Figure 64, the system diverges in the 
opposite direction. In the second case, the system may have been at a higher energy 
level, the system possessing more momentum as it unloaded, overshooting the neutral 
stability point and encouraging divergence in the opposite direction. In both of these 
cases, when the system hits the nonlinear region again, sufficient energy has been added 
to the system so that the hard stops of the torsional spring are hit; oscillations from one 
stop to the other ensue. 

Similar data is presented pertinent to configuration # I . For this configuration, the onset 
of the instability is gentler than that observed for configuration #2. Additionally, the 
system does not hit the hard stops of the spring until after several oscillations from one 
nonlinear range to the other. 
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Figure 67 Nonlinear system divergence for configuration 1 at constant dynamic pressure 

Results for configuration #3 are shown in Figure 68. This are particularly interesting data 
because the nonlinear system is barely unstable and appears to damp slightly instead of 
gaining energy with each oscillation. Prior to 15 seconds, the airfoil is sitting at 1 1 W’ 
angle of attack. The dynamic mode oscillations grow in amplitude, pushing the airfoil 
deeper into the stall region. At 15.75 seconds, the system destabilizes in a static sense. 
The angle of attack changes monotonically until the nonlinear aerodynamic region on the 
negative side is encountered. The nonlinear system diverges again to the positive side 
and then returns again to the negative side. Unlike in the previous cases, the system does 
not have enough momentum to oscillate to the hard stops of the spring. Instead, the 
motion damps slightly. By the second time the system hits the negative side, the 
momentum is no longer sufficient to propel the airfoil through zero. As in the first case 
discussed, it is now equally likely to diverge to either the positive or the negative side. 
Returning to the negative side, the airfoil settles again at the stalled angle of 1 V/ 2 0 . 
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Figure 68 Nonlinear system divergence for configuration 3, tab 976, as velocity is slowly increased 
until divergence condition 


Dynamic mode characteristics at the instability condition 

The time history data presented previously pertained to the system when the rigid angle 
of attack was very close to zero degrees. The characteristics of the system at the 
instability condition appear different when the airfoil is set at a substantial nonzero rigid 
angle of attack. Data is presented in Figure 69 as the configuration #2 destabilizes at 
three different rigid angles of attack. 
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Figure 69 Instability encountered for several values of rigid angle of attack. 
Configuration #2 


The oscillatory behavior appears in the data more prominantly for the non-zero rigid 
body angles of attack, 5° and 6° versus 0°. It is important to note that setting the rigid 
angle of attack at a non-zero angle causes the system to diverge at a lower dynamic 
pressure. Each of the lines plotted in Figure 69 is for a different dynamic pressure. The 
data at 5° angle of attack was acquired at a dynamic pressure of 3. 13 psf. The modal 
frequency is measured as 6.0 Hz. For the same dynamic pressure, the dynamic mode 
frequency of the system at 0° was measured between 5.7 and 6. 1 Hz, depending upon the 
data reduction technique. The frequency of the system at instability for the 6 rigid 
angle of attack is 6.2 Hz. The dynamic pressure is 2.55 psf. At a dynamic pressure of 
2.55 psf, the 0° rigid angle of attack data yielded a frequency between 6.3 and 6.5 Hz. 

The nonzero angle of attack data must be considered in light of the nonlinear regions 
discussed above. There is additional complexity in the system because the structure’s 
equilibrium point is different than the aerodynamic stall equilibrium point. Consider the 
above data for a rigid angle of attack of 5°. The unstable behavior initiates at 
approximately 0.5 seconds as the dynamic mode oscillations grow in magnitude. As time 
progresses, the oscillations are shown to cross the zero angle of attack, without moving 
into the nonlinear effects on the negative side, and then returning to the positive stalled 
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side. This behavior exists due to the disparity in the structural and aerodynamic neutral 
stability points. As the airfoil oscillates towards 5°, the system is unloading structurally 
and aerodynamically. Once it passes 5°. it is becoming structurally loaded such that the 
torsional spring force acts to pull it back towards 5°. The aerodynamic forces are still 
dissipating until the system oscillates beyond 0°. This data set was acquired substantially 
below the instability dynamic pressure of the system at 0° rigid angle of attack. The 
structural restorative moment capability exceeds the aerodynamic moment imparted at 
the low angles of attack. The system is therefore pulled towards the structural 
equilibrium point, 5°. The system oscillates back to the stall angle of attack as a larger 
aerodynamic moment is produced by the airfoil at a now larger angle of attack. 

Similar data is presented for configuration #1 and # 3 in Figure 71 and Figure 70. 
respectively. 



Figure 70 Instability encountered for several values of rigid angle of attack, each at a different 
dynamic pressure, Configuration #1 
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Figure 71 Instability encountered for several values of rigid angle of attack, each at a different 
dynamic pressure, Configuration #3 

Frequency ratios for the wind tunnel model configurations are compared in Table 13. 
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2 
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3 
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0.01 


Table 13 Frequency ratios of experimental data 
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Discussion of Experimental Methods 


Divergence Prediction Methods 

Divergence onset determination using the load monitoring method requires that data be 
acquired very near divergence. Interpretation of data taken well below divergence and 
examined in this way does not easily yield the divergence condition. When using load 
monitoring, it appears best to use a small, nearly zero rigid angle of attack. The load is 
proportional to the elastic increment on angle of attack. This quantity is itself an 
amplification of the rigid angle of attack. The amplification factor becomes very large, 
theoretically goes to infinity, near divergence. As observed in the data for larger rigid 
angles of attack, there is a steady, gradual rise in the moment. This is the effect of 
applying an increasing dynamic pressure to a lifting surface at a non-zero angle of attack. 


Modal Characteristics Determination Methods 

Several methods were utilized to identify modal characteristics. It is not possible to 
completely isolate a single aeroelastic mode’s behavior in an experimental setting. Thus, 
the measurements of the dynamic mode properties contain some content from the static 
mode of the system. 

The logarithmic decrement method is a simple way to calculate damping for systems 
where the damping is low. The time histories produce consistent values in these cases, 
where many cycles can be included inherently in the data processing. The air off data 
presented is a good example of this situation. More highly damped systems, however 
produce a limited number of oscillations. In the case of the plexiglass trailing edge 
configurations, there were often only one or two cycles of decaying motion to analyze, 
even at very low airspeeds. 

It was difficult to extract information using the logarithmic decrement method in the 
vicinity of divergence. Taking several data sets or using different segments of a time 
history produced very different values of damping. The overall system response is a 
combination of all of the system modes. Thus the time history’s damping is not a modal 
damping. Logarithmic decrement results can only be interpreted as modal damping if a 
single mode is dominating the response. Near divergence, it is speculated that two modes 
are contributing- a stable dynamic mode and a barely stable real mode. The calculated 
damping could be the damping of one mode of the other, or more likely could indicate a 
value in between them. 

To utilize the logarithmic decrement technique, a disturbance must be applied to the 
mode. As divergence is neared, the force with which a pluck is administered must be 
reduced. Up to a certain dynamic pressure a very severe pluck can be used. This 
produces more usable cycles of data. In these subcritical data sets, it is speculated that 
the dominant damping effect is produced by the dynamic mode. Very near the 
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divergence condition, it was difficult to pluck the model without destabilizing the system. 
The dominant damping effect here is speculated to be that of the static mode of 
aerodynamic origin which is very near instability. 

The analysis of the logarithmic decrement data taken alone might lead one to consider 
this a case of single degree of freedom flutter, where the dynamic mode destabilizes. The 
additional data and analyses, however, do not support this interpretation. 

An advantage to analyzing the data with frequency domain techniques is that the 
transformation process sorts the information by frequency, separating the information by 
mode. The damping and frequency information pertinent to one mode is distinct from the 
others, as long as the modes are well separated. The main driver in testing configuration 
3 was to raise the frequency of the dynamic mode at divergence so that the peak could be 
fully distinguished. The modal frequency of configurations 1 and 2 was low enough near 
divergence, that it is speculated that there could be substantial “leakage” from the static 
information into the of the dynamic mode measurements. 

The frequency domain techniques encountered difficulties near divergence. At velocities 
well below divergence, a frequency sweep to the gust vanes provided the system with 
sufficient excitation to give clean peaks in the frequency domain. There are several 
issues to overcome in applying these methods near divergence. The first is how to impart 
enough energy into the dynamic mode of the system. The frequency sweeps were 
effective to a dynamic pressure that is very near the soft or linear divergence condition. 
Above this, sine dwell excitations had to be employed. 

In doing the sine dwells, the tunnel condition had to be maintained over along period of 
time. Tunnel drift is thought to have produced a change in dynamic pressure, affecting 
the consistence of a set of sine dwell data. Changes in the dynamic pressure affect not 
only the amplitude of the response at a given frequency, but also shifts the modal 
frequency. Inability to hold tunnel condition is thought to be partly responsible for the 
data scatter associated with analysis of the sine dwell data. 

An additional issue with the frequency domain results is extraction of the important 
parameters from the plot. The faring of the line through the data is subjective and 
becomes more crucial for highly damped systems. 
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CHAPTER FOUR 


COMPARISON OF ANALYSIS AND EXPERIMENT 

The analytical and experimental results, previously presented in this paper, are compared 
for each of the three wind tunnel model configurations. The subcritical eigenvalues of 
the dynamic modes are compared, as well as divergence conditions. 

Configuration # 2, which is described in detail in Table 4, had the lower pitch stiffness 
and the lower pitch inertia. This configuration is considered first. The imaginary part, 
frequency, of the eigenvalue associated with the dynamic mode is plotted as a function of 
dynamic pressure in Figure 72. Analytical results are indicated by the small “x’s”, while 
experimental data are indicated by the larger symbols. Experimental data acquired using 
different methods are represented by the different symbols. Five modes are shown from 
the analysis. Two complex aerodynamic modes are shown; they originate at low 
frequency and spring rapidly to a high frequency. They do not play a role in determining 
either the stability nor the subcritical characteristics of the dynamic mode of interest. 

Two real modes that originate in the aerodynamic model are overplotted on the real axis. 
The mode which originates as the structural dynamic mode starts at the natural frequency 
of the torsion mode, 49.5 rads/sec and migrates to a lower frequency as the aeroelastic 
coupling comes into play. The analytical results and experimental data agree very well. 
Both analysis and experiment indicate that the frequency of the dynamic mode is 
substantially non-zero at the divergence condition. The divergence dynamic pressure 
predicted by the analysis is 4.6 psf. The wind tunnel model diverged at 5. 1 psf and hit 
the hard limit instability at 5.6 psf. 
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Imaginary Part versus versus Dynamic Pressure 11-Mar-2000 



Dynamic Pressure (N/m 3 ) 


Figure 72 Comparison of analytical and experimental results for configuration # 2, frequency of 
system behavior as a function of dynamic pressure 

The real part, indicative of the damping, of the dynamic mode eigenvalue is plotted as a 
function of dynamic pressure in Figure 73. Five modes are shown from the analysis. 
Two of the modes originate at zero and drop rapidly. These are complex aerodynamic 
modes which become very highly damped. The mode which destabilizes is a real mode 
which originated in the aerodynamic model. The mode which originates as the structural 
dynamic mode starts near zero and monotonically decreases throughout the dynamic 
pressure range presented. The fifth mode present from the analysis originated as the 
second real aerodynamic eigenvalue. The subcritical trends in the experimental curves 
follow the trend of the dynamic mode. Near divergence, between 4 and 5.6 psf, there is 
substantial scatter in the experimental determination of damping. The damping values 
generally lie between the analytical value for the dynamic mode and the stable 
aerodynamic mode. While the experiment and analysis are not in perfect quantitative 
agreement, both indicate that the dynamic mode is stable at the divergence dynamic 
pressure. 
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Real Part \ersus Dynamic Pressure 11-Mar-2000 
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Figure 73 Comparison of analytical and experimental results for configuration # 2, damping 
behavior as a function of dynamic pressure 

The real and imaginary parts are combined in a root locus plot in Figure 74. For clarity, 
the complex aerodynamic modes have been removed from the chart. As the reduced 
velocity or dynamic pressure increases, the data progress generally downward and to the 
left. 
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Root Locus as a function of Dynamic Pressure 11-Mar-2000 



Real Part 


Figure 74 Comparison of analytical and experimental results for configuration # 2, root locus as 
dynamic pressure is varied 

The dynamic pressure and frequency results for configuration # 2 are summarized in 
Table 14. 



Table 14 Comparison of analysis and experimental values for configuration # 2 
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Configuration # 1 is considered next. Again, the analytical results and experimental data 
agree well. In the vicinity of the divergence dynamic pressure, the experimental 
frequency increases, departing from the analytically predicted frequency. Both analysis 
and experiment indicate that the frequency of the dynamic mode is substantially non-zero 
at the divergence condition. The divergence dynamic pressure predicted by the analysis 
is 4.6 psf. The wind tunnel model diverged at 5. 1 psf and hit the hard limit instability at 
5.5 psf. 
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Figure 75 Comparison of analytical and experimental results for configuration # 1, 
imaginary part as a function of dynamic pressure 

The real part, indicative of the damping, of the dynamic mode eigenvalue is plotted as a 
function of dynamic pressure in Figure 76. For subcritical conditions below 4 psf the 
shapes of the experimental data agree very well with the analytical dynamic mode. Near 
divergence, between 4 and 5.6 psf, the measured damping values appear to follow the 
real root of aerodynamic origin, which destabilizes. 


Imaginary Part versus versus Dynamic Pressure 10-Mar-2000 
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Real Part \ersus Dynamic Pressure 10-Mar-2000 
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Figure 76 Comparison of analytical and experimental results for configuration #1, real part as a 
function of dynamic pressure 

The real and imaginary parts are combined in a root locus plot in Figure 77. As the 
reduced velocity or dynamic pressure increases, the data progress generally downward 
and to the left. 
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Root Locus as a function of Dynamic Pressure 10-Mar-2000 



Figure 77 Comparison of analytical and experimental results for configuration #1 , root locus as a 
function of dynamic pressure 


The dynamic pressure and frequency results for configuration # 1 are summarized in 
Table 15. 



Dynamic Pressure 

Frequency 


(psf) 

(rads/sec) 

(Hz) 

Analysis: 




Air-off 

characteristics 

0 

21.2 

3.4 

Divergence 

4.6 

6.2 

1.0 

Experiment: 




Air-off 

characteristics 

0 

21.4 

3.4 

Linear System 
Divergence 

5.1 

9.0 

1.4 

Hard Limit 
Instability 

5.5 

13.2 

2.1 

Southwell method 
results 

5.1 




Table 15 Comparison of analysis and experimental values for configuration # 1 
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Configuration # 3 is now examined. The frequency results, shown in Figure 78 
demonstrate the best agreement between theory and experiment of the three 
configurations. The analytical results and experimental data lie on top of one another 
throughout the dynamic pressure range, except right at divergence. In the vicinity of the 
divergence dynamic pressure, there is a small amount of scatter in the measured 
frequencies. The experimental values mainly fall slightly below the analytical 
calculations. Both analysis and experiment indicate that the frequency of the dynamic 
mode is approximately 45 radians/second (7.2 Hz) at the divergence condition. The 
divergence dynamic pressure predicted by the analysis is 14.3 psf. The wind tunnel 
model diverged at 15.2 psf and hit the hard limit instability at 16 psf. 
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Figure 78 Comparison of analytical and experimental results for configuration # 3; frequency versus 
dynamic pressure 

The real part, indicative of the damping, of the dynamic mode eigenvalue is plotted as a 
function of dynamic pressure in Figure 79. Three modes are shown from the analysis. 
The mode which destabilizes is a real mode which originated in the aerodynamic model. 
The mode which originates as the structural dynamic mode starts near zero and 
monotonically decreases throughout the dynamic pressure range presented. The third 
mode shown originated as the second real aerodynamic eigenvalue. The trends in the 
experimental curves follow the trend of the dynamic mode. In the vicinity of divergence, 
between 13 and 14.6 psf, the damping values could be interpreted as being representative 
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of any of the three modes shown from the analysis. Both anlaysis and experiment 
indicate that the dynamic mode is stable at the divergence dynamic pressure. 



0 2 4 6 8 10 12 14 16 18 20 

Dynamic Pressure (psf) 

Figure 79 Comparison of analytical and experimental results for configuration # 3; damping 
characteristic versus dynamic pressure 

The real and imaginary parts are combined in a root locus plot in Figure 80. For clarity, 
the complex aerodynamic modes have been removed from the chart. As the reduced 
velocity or dynamic pressure increases, the data progress generally downward and to the 
left. 
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Figure 80 Comparison of analytical and experimental results for configuration # 3; root locus as 
dynamic pressure is varied 

The dynamic pressure and frequency results for configuration # 3 are summarized in 
Table 16. 



Dynamic Pressure 

Frequency 


(psf) 

(rads/sec) 

(Hz) 

Analysis: 




Air-off 

characteristics 

0 

87.3 

13.9 

Divergence 

14.3 

46.5 

7.4 

Experiment: 




Air-off 

characteristics 

0 

88.0 

14.0 

Linear divergence 

15.2 

45.2 

7.2 

Hard instability 
point 

16.0 

56.5 

9.0 

Southwell method 
results 

15.0 




Table 16 Comparison of analysis and experimental values for configuration # 3 
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CHAPTER FIVE 


SUMMARY OF RESULTS 


Aero elastic Analysis Results 

Discrete time eigenanalyses of the aeroelastic systems revealed a static instability 
originating in an aerodynamic mode and also the characteristics of a noncritical structural 
dynamic mode for the aeroelastic systems examined. The analytical calculations of the 
divergence dynamic pressure agreed exactly with those predicted by the static 
equilibrium equations. 

A database generated by varying the relevant nondimensional parameters revealed that 
the variation of dynamic mode frequency and damping at the divergence condition was a 
complex function of elastic axis position, radius of gyration and mass ratio. The 
dependence on torsional mode frequency was shown to be merely one of scaling. The 
parametric variations revealed that there are regions within the parameter space where 
divergence occurs and the dynamic mode becomes a real mode. The predominant 
category of behavior, however, for the family of configurations studied was shown to be 
a persistence of the dynamic mode, originating from the structural dynamics, at a non- 
zero frequency as the system diverges. 

Additional insight into the modeling and physics associated with system behavior can be 
gained by examining the eigenvectors. The aerodynamic eigenmodes contain the 
essential information of the aerodynamic model. Two static aerodynamic modes exist 
which resemble the static pressure distribution over the wing elements. The complex 
aerodynamic modes are oscillatory wake modes, each at a constant frequency. 

The dynamic mode eigenvector of the aeroelastic analysis yields information regarding 
the stability of the system. The vorticity participation factors can be analyzed to produce 
the modal frequency and damping. The degree of aeroelastic coupling is indicated by the 
amount of vorticity in the wing portion of the vorticity distribution contained in the 
eigenvectors. As velocity increases, the magnitude of the vorticities on the wing 
increases up to a maximum at the divergence condition. 

The static mode of the aeroelastic system that diverges resembles a pressure coefficient 
distribution over the airfoil, similar to the real aerodynamic modes. All of the wake 
vorticity lies in the last element at low velocities. As velocity increases, all of the wake 
elements begin to participate in the mode. Near divergence, the participation of the last 
wake element drops sharply. The sign of the vorticity on this element changes, becoming 
out of phase with the airfoil vorticity and indicating a change in stability of the mode. 



The modal vorticity ratios indicate the relative participation of the aerodynamics and 
structure in each mode. Subcritically, the traditional divergent configurations have 
dynamic modes which contain a large degree of aerodynamic participation. By contrast, 
the configurations studied here that diverge but also have a persistent non-zero frequency 
dynamic mode have modal vorticity ratios that reflect relatively low participation of the 
aerodynamics in that mode. 

Eigenvector orthogonality studies showed that the physical mismatch in eigenvector 
quantities requires that each portion’s orthogonality be considered separately. The wake 
vorticity portion of the dynamic aeroelastic mode eigenvector behave as if the 
aerodynamics are being forced at the modal frequency. Oscillations in the wake 
vorticities show up strongly in the eigenvector phasing if all components of the 
eigenvector are used in the comparison. The vorticity on the wing and the structural 
dynamics considered separately both showed that near zero flow velocity, the dynamic 
mode is nearly orthogonal to the real aerodynamic-originated modes. This reinforces the 
idea that there is no coupling between the structure and the aerodynamics until the system 
is subjected to a substantial flow velocity. This orthogonality is quickly lost as the 
airspeed is increased. As the velocity approaches divergence, the angle between the 
dynamic mode and the unstable mode changes. The modes start to lose their 
orthogonality, allowing energy transference between the modes. 


Experimental Results 

All three experimental configurations diverged in the wind tunnel. The divergence 
dynamic pressure for configurations 1 and 2 was predicted by analysis to be identical and 
was measured to be nearly the same. The physical difference between these two 
configurations was that the first had the Tungsten trailing edge sections, while the second 
had the Plexiglass trailing edge. Classical steady divergence analysis and the analysis 
presented in this study indicate that the divergence dynamic pressure is independent of 
the inertial characteristics of the system. The experimental data agrees with this 
conclusion. Instability of each system was encountered at a slightly higher dynamic 
pressure than the linear theory predicted for divergence. 

The three configurations tested all show the continued presence of the structural dynamic 
originated mode at non-zero frequency when the system becomes statically unstable. 

This was indicated by the time histories as well as frequency domain analysis. 

The three configurations can be compared in terms of the instability onset. The onset of 
instability becomes more sudden and violent for typical sections with low inertias. This 
was experienced in the testing and can be observed by comparing time histories of 
destabilization, as well as by examining the divergence onset prediction results. For the 
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first configuration, with large torsional inertia provided by the Tungsten trailing edge, a 
gentle onset was encountered. 

The destabilizing data at nonzero rigid angles of attack provide compelling evidence that 
the behavior of the dynamic mode is not a strong function of the static stability of the 
system. The frequency of the motion as the system destabilizes is close to the frequency 
that the mode possesses at zero angle of attack for the same dynamic pressure. 

Data was also acquired for test conditions in excess of the linear system divergence 
dynamic pressure. At an angle of attack, between 8° and 1 IV 2 0 , the airfoil reached 
aerodynamic stall. An effect of stalling the airfoil was a reduction of the effective 
aerodynamic moment; dramatic yet understandable nonlinear behavior was thereby 
produced. 



CHAPTER SIX 


CONCLUSIONS & SUGGESTED FUTURE WORK 

The analyses and experiment presented show that aeroelastic divergence can occur 
without a structural dynamic mode losing its oscillatory nature and becoming static. The 
aeroelastic coupling of the static aerodynamic and structural properties that produces 
divergence does not require the dynamic system behavior to cease. Aeroelastic changes 
in the dynamic mode behavior are governed not only by the stiffness, but by damping and 
inertial properties. 

Typical section analysis and wind tunnel experiments demonstrated divergence, 
destabilization in a static sense, but at the same time demonstrated that a dynamic mode 
was still present in the system. These analytical and experimental results challenge the 
basic assumption that divergence occurs as a structural dynamic mode becomes static. It 
has been demonstrated that utilizing dynamic mode tracking to predict divergence onset 
experimentally is inadvisable. The combined aerodynamic and structural stiffness is 
shown to go to zero, but the dynamic mode frequency is shown to not necessarily 
disappear as the divergence condition is reached. 

From this simple analysis and experiment, many possibilities open up for future research. 
In this work, a typical section with a single pitch structural freedom was employed. 
Suggested future investigations include extending the structural degrees of freedom to 
include the plunge mode. Inclusion of the plunge mode could simulate inclusion of rigid 
body or fuselage plunge motion. Study of a wing configuration would be a logical and 
useful extension also. 

The nonlinear behavior which was observed in the present experiment also provides an 
opportunity for follow-on work. A more rigorous investigation of the behavior from a 
phase plane and energy level standpoint might offer interesting results. Additionally, a 
theoretical investigation of the stall behavior is suggested, perhaps utilizing the ONERA 
aerodynamic stall model 1 . 

The analytical method of utilizing a discrete time aerodynamic model could potentially 
be extended to include doublet lattice aerodynamics. As a workhorse in the aeroelastic 
community, an analysis using doublet lattice in this fashion would provide direct 
comparison with current common analysis practices. Insight into the differences between 
roots produced by eigenanalysis and Pade approximation could be gained. 

Experimentally, a simple extension of this work would be to make incremental changes 
in inertial properties in order to produce configurations where traditional divergence 

1 Tang. D. M.. and E. H. Dowell, Comments on the ONERA Stall Aerodynamic Model and its Impact on 
Aeroelastic Stability 
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occurs. The current model would require relocation of the axis of rotation closer to the 
center of pressure and addition of mass, external to the airflow. 

Acquisition and processing of unsteady pressure data offers additional research 
opportunities. The analytically determined eigenvector study has offered some insight 
into the phase relationship between the structural displacement and velocity and the 
vorticity distribution. The vorticity or pressure distribution, measured in an experiment, 
could provide additional insights. Potentially, pressure data could be examined and 
proper orthogonal decomposition techniques applied to examine aerodynamic modal 
participation in overall system response. 
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APPENDIX A 
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EQUATIONS OF MOTION 
NOMENCLATURE 

Aerodynamic matrix 

Distance from midchord to elastic axis 

Aerodynamic matrix 

Semi-chord 

Aeroelastic force matrix 

Intermediate calculation matrix, defined in equation 46 
Intermediate calculation vector, defined in equation 47 
Pitch mode damping 
Structural dynamic matrix 
Downwash matrix 

Elastic axis position, measured positive aft from the center of 
pressure 

Aerodynamic load vector 
Integral expression, see equation 40 
Torsional inertia 

Intermediate structural dynamic matrix, defined in equation 54 
Intermediate structural dynamic matrix, defined in equation 54 
Number of aerodynamic elements on the wing 
Mass 

Total number of aerodynamic elements 
Time step number 
Generalized structural coordinates 
Dynamic pressure 
Radius of gyration 

Intermediate calculation vector, defined in equation 44 
Intermediate calculation vector, defined in equation 45 
Velocity 

Reduced velocity, (V=U/co H b) 

Downwash 

Vector of locations of vortices in the aerodynamic model; chord- 

wise location 

Discrete time eigenvalue 

Time step size, temporal discretization 

Aerodynamic element size, spatial discretization 

Vorticity vector 

Aerodynamic kernel function 

Torsional stiffness 
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M Moment 

(X Aerodynamic relaxation factor 

a Angle of attack 

\ Continuous time eigenvalue 

p Mass ratio 

p Density of air 

( 0 U Frequency of torsional mode (radians/second) 

q Vector of locations of collocation points in the aerodynamic model 

£ Damping 

Superscripts 

n time step number 

Subscripts 

a Pertaining to the torsional degree of freedom 

oo Freestream 

0 Steady quantity 

1 Unsteady quantity 

a Airfoil motion quantity 

EA Quantity at the elastic axis 

i,j Designation for an element of a matrix which lies in the ith row, 

jth column 

k kth aerodynamic element 

M The Mth aerodynamic element 

wake Quantity in the wake 

wing Quantity on the wing 


Aerodynamic equations 

The aerodynamic model was constructed by considering time-marching relationships as 
the vorticity develops, on the wing and in the wake. For this study, the airfoil is modeled 
as a 2-dimensional flat plate. The airfoil and the wake are divided into segments, referred 
to as aerodynamic elements. 

The vortex lattice aerodynamics are generated by placing vortices of strengths to be 
determined at points on the airfoil and in the wake. Control points, usually located aft of 
the vortex locations, are points where the boundary conditions must be satisfied. Typical 
placement is for the vortices to be located at the '/4-chord point of an aerodynamic 
element. The control points are typically placed at the 3 /4-chord locations of the elements. 
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There are 3 basic relationships contained in the vortex lattice equations, detailed in the 
following paragraphs, which are combined to form a matrix equation. 

The velocity induced by the discrete vortices is set equal to the downwash caused by the 
airfoil’s motion. Equation 24. 

•ff ' = Equation 24 

j = I 


i = 1 M 

The kernel function, Kij, relates the vorticity at point j, Tj, to the downwash generated at 
point i, Wj. For an isolated flat-plate airfoil in two-dimensional incompressible flow is 
given in Equation 25. 

1 

El) = — r Equation 25 

2n\xi -g) 

Applying Kelvin’s theorem generates the second basic relationship. Quoting Hall 1 
“unsteady vorticity is shed into the wake; its strength is proportional to the time rate of 
change of circulation about the airfoil. If the time step is taken to be equal to the time it 
takes the vorticity to convect from one vortex station to the next, then the strength of the 
first vortex point in the wake at the time n+1 is given by (Equation 26)” 

+ l = — X (r/ ,+ ' — Ty ) Equation 26 

7 = 1 

Once the vorticity has been shed into the wake, it is converted in the wake at the 
freestream velocity. Convection provides the final relationship utilized in constructing 
the aerodynamic equations. Fixing the time step such that At = U Ax, this convection is 
described by Equation 27. 

= P p j Equation 27 


,i = (M +2),(N - 1) 


Because the wake is modeled with a finite length, the convection relationship for last 
vortex element must be treated specifically. “Otherwise, the starting vortex would 
disappear abruptly when it reached the end of the computational wake, producing a 


1 Hall, Kenneth C., Eigenanalysis of Unsteady Flows about Airfoils, Cascades and Wing 
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discontinuous change in the induced wash at the airfoil. To alleviate this difficulty, ... the 
vorticity is allowed to dissipate smoothly by using a relation factor, (Equation 28).” 


T-i n + 1 r' /i 

1 i ~ 1 /-I 


+ a r? 


Equation 28 


,i = N 


The equations are combined into the matrix expression. Equation 29. The Kernel 
function relationship between the vorticities on the wing and wake and downwashes on 
the wing form the first M rows of the equations. Kelvin’s theorem is seen as the (M+l ) 
row. Convection in the wake appears in rows M+l through N. 
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-a 
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Equation 29 

Rewriting Equation 29 in general terms produces Equation 30. 
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Mrr 1 + w =vr 


Equation 30 


Aerodynamic moment equations 

The moment about the elastic axis generated by the aerodynamics is calculated using the 
integral expression in Equation 31. The distribution of vortex strength per unit length in 
the chord-wise direction is specified as y(^). The distance from the elastic axis to the 
mid-chord of the typical section, measured negative aft, is specified as ab. 

Mg A (t) = p(x- ab\uy{x,t) + ^\* b y(^,t)d^\^X Equation 31 

The moment equation is spatially discretized; approximations to the spatial integrals are 
made. In doing so, the vorticity, y(^k), of each aerodynamic element is used. The 
moment arm for the force generated by each element’s vorticity is the distance from the 
vortex location to the elastic axis. 

The integral is broken apart into the steady and unsteady portions, denoted M 0 and Mi, 
respectively as shown in Equation 32 and Equation 33. 

Mq ( t ) = \ h _ h p(x - ab py(X,t )dx Equation 32 


M ] {t) = \ h _ b p{x-ab 

The moment at a given time, (n+l/2)At, is the sum of the steady and unsteady parts at that 
time. Equation 34. 

M ea (( w + / K)a/)= M o ((« + M \ ((/I + Equation 34 

It is desired to express the moment in terms of the vortex strength distribution at integer 
time steps. Before proceeding, the notational convention exhibited in Equation 35 is 
introduced. The superscript represents the time step and the subscript represents the 
spatial element location. 

yig k , (n + y 2 )\t ) = y n k + ^ 2 Equation 55 

The vortex strength at a location on the wing is then approximated at time (n + Vi )At as 
the average of the value for the preceding and following time steps. Equation 36. 


j t \ 'IbftZ’MZ dx 


Equation 33 
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Equation 36 



H-ry+rQrr' 

? 


The steady portion of the moment is approximated using a weighted summation of the 
vortex strength distribution on the airfoil as shown in Equation 37. The vortex 
associated with the k' h aerodynamic element is located at 4 . This is the location relative 
to the midchord and ab is the distance from the midchord to the elastic axis. Note that the 
control or downwash point of each aerodynamic element, x k , is located half an element 
aft of the vortex location, x^= 4 + V 2 Ax. The moment arm for each vortex utilizes the 
vortex location. The moment computation, however, is evaluated for the downwash 
locations. This seeming incompatibility is necessitated by the equations to which the 
moment expression is coupled in forming the governing aeroelastic system equations. 


1/ M 

M () " + /2 =-UpAx 1(4 -ab)Yk 

k = \ 


Equation ,? 7 


The summation is implemented as a vector product as shown in Equation 38. 


Mq —U p&x\(£\ — (tb) (£2 ab) (4/ \l [^ wing I* 


Equation 38 


The unsteady portion of the moment equation, Equation 39, is next evaluated at the time 
under consideration. In addition to the vortex strengths varying with time, there is a 
derivative expression which must be considered. 


M 



\ b _ b p{x-ab 




-1 n+\ 


dx 


Equation 39 


The internal integral at time n, evaluated at the downwash point of the k th aerodynamic 
element is represented as I|, n . The integral takes on a different value for each 
aerodynamic element, as the upper limit of integration is the chord-wise location of the 
control point for a given element. The integral expression can be approximated as a 
summation of contributions from each aerodynamic element, as shown in Equation 40. 
The contribution from each element is the elemental vortex strength multiplied by the 
length of the element which is upstream of the downwash point under consideration. 
That is, if an element is ahead of the element under consideration, its entire length is 
utilized in the calculation. If an element is behind the element under consideration, none 
of its strength is used in the calculation, so it is multiplied by zero. Three-fourths of the 
element being considered lies upstream of the downwash point, so % of that element's 
length is utilized in the summation. The summation is also presented in the form of a 
matrix multiplication suitable for implementation. 
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= n y{ £ )d z = (x ^ i y] = jt-o 

7 = 1 
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/4 





Equation 40 

The time derivative of the integral expression is approximated using the central 
difference theorem, Equation 4 1 . 


(dTi 

dt 


n+ 


r A/_ 


n + 


r rt + 1 




A/ 


Equation 41 


The unsteady portion of the moment is approximated by a double summation, Equation 
42, where one of the summations is included in the evaluation of the integral expression. 


M 


n+ 


— pArf '1 


* = l 


Equation 42 


This can be written as a matrix product, Equation 43. 


M" ^=pU^-ab) is2 -ab) - (§ M -ab ) j 


/ 4 

1 ) 


0 ••• 0 
0 

i x 


tty. + >) 

VI/ wing J V wing J / 


Equation 43 


Implementing the calculations for the steady and unsteady moments can be combined 
utilizing the following matrix products. Equation 44 through Equation 47. 


*l=yLC£|-«*) (f 2 ~ ab ) ■■ (Zm -ab)\ 


Equation 44 


h = pL(£i ~ ab ) (£2~ab) 


3/ 
/ 4 


(£m ~ ab \ ■ 


o 

3/ 

/4 



Equation 45 


Zero padding is required to make the dimensions such that the matrices fit into the 
aeroelastic system equations. 
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Equation 46 


c i 


goo °(i,^} °(i,jv) 

mb 2 bl-h) °(\,N-M) 


c 2 =- 



%N e ) 

(*1 + * 2 ) 


°{1JV-A/) 


Equation 47 


The dimensions of the coupling equations necessitate the formation of a vector 
expression where the aerodynamic moment equation is the second row. 


/ 


n+/ 2 


j 0 l" + ^ 

\Mea\ 


Equation 48 


A vector containing all of the vortex strengths is defined in Equation 49. 


r" = J 

1 f. 

1 wing _ 


Vr 

> = < 

7\ 

w 

\v n 

l wake . 


f N 


JN . 


Equation 49 


The general form of the aerodynamic moment equation is then shown in Equation 50. It 
is noted again that the first row is a zero element and only the wing elements of the 
vorticity distribution are utilized in computing the aerodynamic moment. 


y ,,+ >2 =c 2 r” +1 +c,r” 


Equation 50 


Structural Dynamic Equations 

The equations of motion for a typical section with a single pitch degree of freedom, 
possessing inertia, damping and stiffness characteristics, and subjected to an aerodynamic 
moment are given in Equation 51. 

I C K M y Equation 51 

To write the equation in terms of nondimensional quantities, Equation 5 1 is divided by 
mass and semi-chord squared. 
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Equation 52 


r a fc}+ W 2 H+ r a <°l M= My / . 2 

/mb /mb 

Equation 52 can be rewritten as a set of first order equations. 


'1 0 ■ 

[ a) 


< > + 
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0 -1 

r 2 co 2 ^ a > 
r a UJ a / , . 

mb 


a\ 

a\ 


M/ -1°) 


The above matrices and vectors are redefined in Equation 54. 


Equation 53 


if 2 } + \J I 1+ {f } Equation 54 

The time derivative of the generalized structural coordinate at a given time, (n+l/2)At, 
can be expressed approximately using the central difference theorem. Equation 55. 


fer'-w 

dt\ 1 At j A t 


Equation 55 


The generalized structural coordinate can be approximately using the central difference 
theorem. 

= y 2 (fef +l + -fcrf ) Equation 56 

The first order equations can be written for a specific instant in time. Writing them at 
time (n+ Vt )At and employing the central difference approximations shown in Equations 
55 and 56 leads to Equation 57. 


l2 + h~ 

fer' + 

\~-h + J .1 

CN 


At 2 _ 


w+{/r>M>} 


Equation 57 


This is now in the general form represented in the main text of the paper. 
D 2 q n+ '+Diq n +f n+l = 0 


Equation 58 


Downwash Equations 


The downwash is defined as the vertical component of velocity on the airfoil, w a , positive 
downward. Using the velocity potential function, <|>, the downwash can be expressed in 
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terms of temporal and spatial derivatives. The vertical displacement, z, is defined as 
positive upward. 


W a ~ ^ Equation 59 

dz dt dx 

The vertical motion of the airfoil, z a , is a function of both spatial coordinate, x, and time, 
t. The motion can be represented in terms of a spatial mode shape, <J>(x), and a temporal 
generalized coordinate, £(t). 

Z a (.* , t ) = 0 (* )£ (0 Equation 60 

The downwash expression can then be written as Equation 61, where the negative signs 
have been incorporated into the mode shapes. 

W a = ®(x)'^p- + U 00 Equation 61 

a dt dx 

The generalized coordinates and their derivatives with respect to time can be written as a 
single vector. 


M= 


#(o 1 

jg(0 

dt I 


The downwash equation can be written in tirst order form. 


Equation 62 


W 


a 


d<S>(x) 

dx 


O(x) 




Equation 63 


Specifically, for the pitch degree of freedom, the vertical displacement of a typical 
section is given in Equation 64. The displacement at a point on the airfoil is a function of 
its location, x, measured relative to the center of rotation. The sign convention utilized is 
positive distance is forward of the center of rotation and angle of attack is positive nose- 
up. 


z a = -xa 


Equation 64 


<t>(x) = — X 


Equation 65 
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£(t ) — Ot Equation 66 

The downwash equation, written for a single point on the wing, in general form is given 
in Equation 67. 

W a =[~Uoo -jrjfer} Equation 67 

Formally, the generalized coordinate vector is comprised of the structural generalized 
coordinate and its first derivative with respect to time. For the case of the pitch only 
typical section the generalized coordinate is the angle of attack. 



Equation 68 


The downwash equations can be written in the general form presented in the main text. 


w n = Eq n 


Equation 69 



APPENDIX B 


EIGENVECTOR INVARIANCE UNDER TRANSFORMATION 

NOMENCLATURE 

A State matrix 

B A general matrix 

C Jordan form matrix 

F Matrix which is similar to B 

P A general matrix 

R Transformation matrix 

s Laplace variable 

t Time 

T Time step size 

X(s) Laplace domain representation of state variable vector 

x(t) Vector of state variables 

y Substitution variable for state vector, see equation 89 

z Unit delay operator, discrete time eigenvalue 

O State transition matrix 

(3 Eigenvalues of B matrix 

X Continuous time eigenvalue 

^ Eigenvector 

SUBCRIPTS 

A, B, C, F Pertaining to matrix A, B, C, or F 
c Continuous time quantity 

d Discrete time quantity 

e AT Pertaining to discrete time matrix, e AT 


Theorem: The eigenvectors for the associated continuous and discrete time state space 
systems are equal. 

Proof. 


• Let the continuous time state space equations for the unforced system be given by 
Equation 70. 


jf(/)= Ax(t ) 


Equation 70 



A solution can be obtained by utilizing the Laplace transform. Equation 71, where X(0) is 
the state of the system at the initial condition. 

sX(s) — ATO) = ALfs) Equation 71 

The system of equations can be rearranged as in Equation 72. 


X{s) = [ls-A]~ ] X(0) 

Define the state transition matrix at a time t,O c (t), as in Equation 73. 


Equation 72 


<f> c (t) = L 




1 ? 1 

I + At H — At H — At 3 + • • 


2! 3! 

Rewriting in summation notation produces Equation 74. 


Equation 73 


<M')= I 


A k t k 


-n k\ 


= e 


. At 


Equation 74 


k= 0 


x(t) = O c ( t )x( 0) Equation 75 

A non-zero initial time can be accommodated in these equations. 


X(t)- <t» c (t — t() )X(fg) 


Equation 76 


-/o) = 


y A k {t-tp) k 

i o k\ 


Equation 77 


• To obtain the discrete time model with a time step size of T seconds, Equation 77 
can be evaluated at time t=nT + T, with to = nT. 


x(nT + T) = <t> c (nT +T — nT)x(nT) 


Equation 78 
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x((n + 1)7) = O c {T)x(nT) 


Equation 79 


In discrete time notation, the time step size is generally omitted. 

x[n + 1 ] = <t> r (T )x[// ] Equation SO 

The discrete time state space equations for the unforced system are given by Equation 81 . 

*(" + 0 = A c j X (n ) Equation S 1 

Comparing Equation 80 Equation 81, the continuous and discrete time system matrices 
are related as in Equation 82. 

oo y A 

A d -e AT = <J> C (T )= X A k Equations 2 

k=0 k ' 

The discrete time equations can thus be written as Equation 83. 

X (/? + 1 ) = e A T X (n ) Equation S3 

• For the continuous time relationship, assume the time variation to be simple 
harmonic motion, Equation 84. 

JC (t ) = e M X (0 ) = e M £4 Equation S4 

An eigenvalue problem is formulated by substitution of the simple harmonic solution into 
the state space equations. Equation 85. 

Equation 85 

• For the discrete time relationships, approximate the time variation as a unit delay. 
Equation 86. The response of the system as an initial time is given a new notation. 

X(n) — z”jr(0) =z n ^ g AT Equation 86 

An eigenvalue problem is formulated by substitution of the unit delay solution into the 
discrete time state space equations. Equation 87. 

AT n 

€ C at = Zs AT Equation 87 

J e e 

• Thus the theorem can be restated as the eigenvectors of A, denoted £, A , are equal 
to the eigenvectors of Aj, denoted ^ e AT . 
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Following the development of Kincaid and Cheney 1 , “if we possess the Jordan form C of 
A and we know the transformation (given in Equation 88), then we can change variables 
by substitution as in (Equation 89).” 

£4 1 A £4 = C Equation 88 


X C = 4 Ay Equation 89 

The differential equation. Equation 70, and the prescribed initial condition, x(0), can be 
recast as Equation 90 and Equation 91. 

q A y = A<; A y Equation 90 

g^y(())= x(0) Equation 91 

Rearranging Equation 90 produces Equation 92. 


y=g A A<^ A y Equation 92 

Substituting the Jordan form matrix, C, from Equation 88 produces Equation 93. 

y = Cy Equation 93 

This set of equations can then be solved, Equation 94. 


y(t )= e CT y(0) Equation 94 

Substituting this solution into Equation 89 and reverting to the original state vector 
produces Equation 95. 

*(/) = % A e CT £; A x(o) Equation 95 

Comparing Equation 95 with Equation 75 leads to the expression of the matrix 
exponential in terms of the Jordan form matrix exponential, Equation 96. 

e A T = e CT £4 1 Equation 96 


1 Kincaid. David, and Ward Cheney, Numerical Analysis: Mathematics of Scientific Computing, page 562 
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The eigenvectors of a general matrix, B, can be obtained from the eigenvectors of any 
similar matrix, F, if the transformation between the two matrices is known. The 
transformation relationship is defined. Equation 97 by the matrix R. 

PR — B Equation 97 

The eigenvalues of a matrix are invariant with coordinate transformation. The 
eigenvalues of both matrices F and B are denoted P in this development. The 

eigenvectors are denoted The eigenvalue problem for B is given by Equation 98. 

Bt, ft =£bP Equation 98 

Substituting Equation 97 and premultiplying by R results in Equation 100. 

R FRg B =£ B /3 Equation 99 

Examining this equation, it is recognized that this is an, eigenvalue equation for the matrix 
F, where the eigenvectors are shown to be R^ B - 

FiR^B )=(/?£ b )/5 Equation 100 


£p =R^g Equation 101 

The preceding derivation is now applied to Equation 96, making the following 
substitutions. 


r -Za 

Equation 102 

o 

II 

Equation 103 

11 

Equation 104 


This results in an expression for the discrete time system matrix in terms of the Jordan 
form system matrix and the continuous time eigenvectors, Equation 105, and a 
relationship among the eigenvectors, Equation 106. 

e AT _ ^£—1 V 1 gCT £-1 Equation 105 
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\CT ~^ A ] ^ e AT 


Equation 106 


C is the Jordan form matrix associated with A, the continuous time system matrix. The 
formula for taking the exponential of a general matrix, P is given in Equation 107. 


PT 

e 



Equation 107 


If P is a matrix of diagonal elements, the P k is a matrix of each element to the k th _jpower. 
Thus for P diagonal, e PT is diagonal and of the same matrix dimension. Thus, e CT is a 
Jordan form matrix. The case of the non-diagonal Jordan matrix has been left to the truly 
inspired reader. 

The eigenvectors of a Jordan form matrix form an identity matrix. 


4 


e 


CT 


= / 


Equation 108 


Using this relationship and Equation 106 leads to Equation 109. 

I =4 A l 4 e AT 


Equation 109 


Therefore, the eigenvectors associated with the continuous and discrete time state space 
equations are equal, Equation 1 10. 

<^4 = £ 1 7 Equation 110 
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